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INTRODUCTION 


An  important  problem  which  must  be  faced  by  the  designer  of  advanced 
gas  turbine  engines  is  the  prediction  of  the  viscous  flow  field  in  and  around 
turbine  and  compressor  blade  passages.  Such  a prediction  is  required  to 
determine  heat  transfer  rates  and  aerodynamic  losses  both  of  which  may  be 
critical  to  successful  engine  operation.  Inaccurate  predictions  of  either 
loss  coefficients  or  heat  transfer  rates  may  result  in  poor  estimates  of 
engine  performance  or  even  catastrophic  failure  of  the  engine  components. 

For  example,  excess  aerodynamic  losses  associated  with  viscous  phenomena 
including  boundary  layer  separations,  trailing  edge-wake  interactions  and 
massive  stall  may  result  in  a serious  deterioration  of  component  efficiency. 
In  addition,  excessive  heat  transfer  rates  associated  with  boundary  layer 
separation  and  reattachment  on  turbine  blades  and  end  walls  can  have  damag- 
ing effects  as  the  resulting  hot  spots  may  result  in  structural  failure. 

Since  aerodynamic  losses  and  heat  transfer  rates  are  associated  with  the 
viscous  nature  of  the  fluid,  the  ability  to  predict  the  viscous  flow  field 
in  high  performance  turbine  and  compressor  blade  passages  i6  crucial  to  the 
successful  design  process. 

For  the  flow  regimes  and  configurations  of  practical  interest,  in  either 
internal  or  external  aerodynamics,  the  Reynolds  number  is  usually  suffi- 
ciently high  so  that  viscous  effects  are  concentrated  in  relatively  thin 
layers.  While  this  in  general  allows  the  direct  use  of  the  boundary  layer 
approximation  to  obtain  meaningful  predictions  of  viscous  phenomena  there  are 
several  regions  of  flow  over  airfoils  where  the  boundary  layer  approach 
breaks  down.  For  example,  this  approach  does  not  apply  in  the  vicinity  of 
leading  edge  separations,  shock  wave/boundary  layer  interactions,  strong 
blowing  sites  and  trailing  edge-wake  interactions.  The  common  feature  of 
the  flow  in  these  regions  is  that  the  boundary  layer-like  viscous  region  is 
displaced  from  the  airfoil  surface  and  exerts  a significant  influence  on  the 
inviscid  flow.  Therefore  it  is  not  surprising  that  the  flow  in  such  regions 
cannot  be  adequately  described  by  the  usual  boundary  layer  approximation. 

When  viscous  displacement  effects  do  alter  the  inviscid  flow  field  signifi- 
cantly, it  is  necessary  for  the  calculation  procedure  to  recognize  the  mutual 
dependence  between  the  viscous  and  inviscid  regions  either  by  a solution  of 
the  full  equations  of  viscous  fluid  motion  throughout  the  entire  region  of 
interest  or  by  a strong  interaction  analysis  between  a viscous  region  solu- 
tion and  an  inviscid  outer  field  solution.  Thus  two  methods  of  approach  are 
currently  available  which  offer  the  long  term  prospect  of  providing  useful 
viscous  aerodynamic  design  information.  One  consists  of  obtaining  numerical 
solutions  to  the  compressible  Navier-Stokes  equations  and  the  other  consists 
of  constructing  and  numerically  solving  viscid-inviscid  interaction  models  of 
localized  flow  regions. 


The  emergence  of  numerical  solutions  to  the  Navier-Stokes  equations  as 
a viable  method  for  predicting  viscous  flows  is  a fairly  recent  phenomena 
resulting  from  rapid  advances  in  numerical  analysis  coupled  with  development 
of  high  speed  computers.  Such  solutions  can  be  obtained  either  through  solu- 
tion of  the  steady-state  equations  via  a relaxaation  procedure  or  through 
solution  of  the  time-dependent  Navier-Stokes  equations  under  the  influence 
of  steady  boundary  conditions.  In  the  latter  case,  each  time-step  may  be 
thought  of  as  a step  in  the  relaxation  procedure.  Time  dependent  solution 
procedures  may  be  either  explicit  or  implicit.  If  the  method  is  explicit, 
the  time-step  is  governed  by  stringent  stability  limits  which  relate  the 
maximum  time-step  to  the  size  of  the  computational  grid.  If  a viscous  layer 
is  to  be  adequately  defined,  a fine  grid  is  required  in  the  vicinity  of  the 
blade  surface  and  in  such  cases  the  stability  limit  would  make  an  explicit 
calculation  impractical.  However,  implicit  methods  are  not  subject  to  such 
stability  limits,  rather  they  are  only  limited  by  the  physical  time  scale  of 
the  flow.  As  a result  a time-dependent  implicit  solution  procedure  for  the 
three-dimensional  Navier-Stokes  equations  has  been  developed  at  UTRC  by 
Briley  and  McDonald.  This  computer  code  is  a highly  modularized  program 
which  has  options  for  two-  and  three-dimensional  modes  of  operation.  In 
two-dimensions,  solutions  to  a general  coordinate  representation  of  the 
Navier-Stokes  equations  can  be  obtained.  The  code  has  been  used  previously 
to  calculate  laminar  and  turbulent  flows  in  ducts  and  low  to  moderate 
Reynolds  number  flows  past  simple  aerodynamic  shapes,  and  it  could  form  a 
basis  for  the  calculation  of  viscous  flows  through  cascades.  Such  solu- 
tions would  include  viscous  effects  and  would  be  extendable  to  turbulent 
flows  and  to  three-dimensional  and  time-dependent  problems.  It  is  antici- 
pated, however,  that  the  computing  time  required  to  obtain  Navier-Stokes 
solutions  for  cascade  flows  may  limit  the  usefulness  of  such  solutions  to 
the  turbomachinery  designer.  In  addition,  there  is  some  doubt  at  present 
that  stable  Navier-Stokes  solutions  can  be  obtained  in  the  Reynolds  number 
range  of  practical  interest  (i.e.,  Re  ■ 0(10^).  Future  research  should 
do  much  to  alleviate  these  limitations  but  at  this  juncture  the  prospect  of 
effectively  using  Navier-Stokes  calculations  for  design  applications  still 
seem  to  be  far  off. 


Although  the  Navier-Stokes  equations  contain  all  the  necessary  physics 
to  describe  viscous  separation  phenomena,  it  has  long  been  felt  that  an 
intermediate  theoretical  approach,  i.e.,  a viscous/inviscid  interaction 
model  could  provide  a useful  basis  for  describing  a large  class  of  boundary 
layer  departure  flows.  The  results  of  relatively  recent  numerical  experi- 
ments do  provide  empirical  support  for  the  use  of  such  models  for  both  low 
and  high  speed  flows.  In  those  cases  where  detailed  investigations  have  been 
conducted,  predictions  based  on  viscous/inviscid  interaction  models  have  been 
found  to  reproduce  Navier-Stokes  predictions  and/or  experimental  data  at 
moderate  to  high  Reynolds  numbers. 
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Viscid-inviscid  interaction  models  (interaction  boundary  layer,  triple 
deck,  etc.)  consider  three  basic  elements  of  the  flow:  an  inviscid  region, 
a viscous  region  which  is  shear  layer  like  and  an  interaction  region  which 
actively  connects  the  two.  Such  a model  allows  one  to  focus  attention  on 
the  significant  phenomena  in  their  respective  flow  regions  and  therefore 
develop  efficient  and  reliable  solution  techniques  tailored  to  the  needs  of 
each  individual  region.  Numerous  examples  exist  where  the  interacting  bound- 
ary layer  concept  has  been  successfully  applied  to  laminar  and  turbulent 
separation  bubbles  in  external  aerodynamics  problems.  In  addition,  the 
triple  deck  equations,  which  are  a subset  of  the  interacting  boundary  layer 
equations,  have  been  used  to  rigorously  study  and  solve  trailing  edge  flows. 
Much  fundamental  work  has  been  done  for  laminar  flow  i.ast  sharp  and  blunt 
(unstalled)  trailing  edges  and  Melnik,  Chow  and  Mead  have  achieved  a dra- 
matic success  for  at  least  one  highly  important  practical  problem,  the  pre- 
diction of  Reynolds  number  effects  on  the  lift  (Kutta  condition)  for  turbu- 
lent flow  past  transonic  airfoils  with  sharp,  attached  flow,  trailing  edges. 
Finally,  recent  fundamental  work  indicates  that  the  combination  of  classical 
free  streamline  theory  and  viscid-inviscid  interaction  concepts  will  provide 
a valid  model  for  predicting  the  major  features  of  massive  stall  phenomena. 

It  currently  appears  that  viscid-inviscid  interaction  models  have  the 
potential  of  providing  fast  and  accurate  predictions  of  several  of  the  vis- 
cous phenomena  of  importance  in  turbomachinery  applications.  As  such,  inter- 
action concepts  appear  to  hold  the  promise  of  producing  useful  and  efficient 
design  calculation  procedures  for  turbomachinery  appl icat ions . Since  the 
interaction  equations  are  a subset  of  the  Navier-Stokes  equations  and  as 
such  identify  the  dominant  terms  in  those  equations  and  the  appropriate 
scaling  laws  and  correlation  parameters,  then  development  of  numerical  tech- 
niques for  viscid-inviscid  interaction  models  may  have  a direct  impact  on, 
and  indeed  may  be  a prerequisite  for  the  development  of  reliable  and  effi- 
cient ways  to  solve  the  Navier-Stokes  equations  in  the  high  Reynolds  number 
regime . 

With  the  foregoing  considerations  in  mind  a research  effort  has  been 
conducted  to  develop  an  analysis  for  the  prediction  of  high  Reynolds  number 
flow  in  a cascade  passage.  Under  the  present  study  an  existing  implicit 
time-marching  Navier-Stokes  computer  code  has  been  further  developed  to  treat 
cascade  flows.  In  particular,  the  general  coefficient,  matrix  inversion,  and 
boundary  condition  subroutines  have  been  modified  to  allow  for  the  specifica- 
tion of  cascade  blade-to-blade  periodicity  conditions.  In  addition,  calcu- 
lations have  been  performed  to  assess  the  reliability  and  accuracy  of  the 
Navier-Stokes  computer  code  and  to  verify  the  successful  completion  of  the 
coding  changes  made  under  the  present  contract.  Finally,  an  effort  ha6  been 
initiated  towards  bringing  viscous/inviscid  interaction  concepts  to  bear  on 
the  viscous  phenomena  of  importance  in  turbomachinery  flows.  In  particular, 
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a viscous/inviscid  interaction  model  incorporating  asymptotic  triple  deck 
concepts  has  been  formulated  for  high  Reynolds  number  laminar  trailing  edge 
flows.  Numerical  solution  procedures  based  on  the  triple  deck  version  of 
this  model  have  been  developed  and  results  for  both  attached  and  separated 
"stalled"  trailing  edge  flows  have  been  obtained. 
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PART  I 

DEVELOPMENT  OF  A NAVIER- STOKES  SOLUTION  PROCEDURE 
FOR  VISCOUS  CASCADE  FLOWS 
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LIST  OF  SYMBOLS  - PART  I 


Except  where  dimensions  are  specified,  all  quantities  in  the  following 

are  nondimens ional ; physical  velocities  are  normalized  by  U , density  by  p_, 

2 r r 
pressure  by  PrUr  , dynamic  viscosity  by  pr,  and  time  by  Lr/Ur 

where  Lf  is  the  reference  length. 

a1-1  • i,  Coefficient  matrix,  Eq.  (1-21) 

A > J > * 

A Constant  defined  in  Eq . (1-3) 

Coefficient  matrix,  Eq.  (I-15a) 

b?  • v Coefficient  matrix,  Eq.  (1-21) 

B Constant,  defined  in  Eq.  (1-3) 

Coefficient  matrix,  Eq.  (1-21) 

Column  vector,  Eq.  (1-21) 

Finite  difference  operators  for  coordinate  ym 
Momentum  equation  coefficient,  Eq.  (1-6) 

Damping  constants  for  coordinate  transformations,  Eqs.  (1-24) 

Column  vector  function  of  the  dependent  variables  and  their 
spatial  derivatives  in  a single  coordinate  direction,  Eq.  (1-9) 

Boundary  condition  matrix  operator,  Eq.  (1-30) 

Column  vector  containing  spatial  derivatives  of  the  dependent 
variables  in  the  ym  - direction  only,  Eq.  (1-17) 

Coordinate  base  vector,  Eqs.  (1-35) 

Momentum  equation  coefficient,  Eq.  (1-7) 

Momentum  equation  coefficient,  Eq.  (1-6) 

Metric  tensor  coefficient 

Inverse  metric  tensor  coefficient 
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LIST  OF  SYMBOLS  - PART  I (cont'd) 

^ijk  Momentum  equation  coefficient,  Eq.  (1-6) 

H Column  vector  function  of  the  dependent  variables,  Eqs.  (1-9,  1-10) 

I Boundary  point  index  - y*  - direction 

Ijj  Momentum  equation  coefficient,  Eq . (1-6) 

J Jacobian 

2 

Boundary  point  index,  y - direction 
K Momentum  equation  coefficient,  Eq.  (1-6) 

Lf  Reference  length,  meters 

L**j  Momentum  equation  coefficient,  Eq.  (1-6) 

J Matrix  operator,  Eq.  ( I — 1 5c ) 

Matrix  operator  containing  spatial  derivatives  in  the  ym  - 
direction  only,  Eq . (1-17) 

Mr  Reference  Mach  number 

n Unit  outward  normal  vector 

p Pressure 

r Computational  radial  coordinate 

R Physical  radial  coordinate 

Rj  Momentum  equation  coefficient,  Eq.  (1-6) 

Re  Reference  Reynoldb  number,  PrUrWur 

Re^x®  Mesh  Reynolds  number,  |UB|Ax*Re 

S Column  vector  containing  mixed  second  order  spatial  derivative 

terms,  Eq . (1-9) 

t Time 
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LIST  OF  SYMBOLS  - PART  I (cont’d) 


T 

T„ 


v(i) 


a 


® 


Y 

fij 

1 

At 


Ax™ 

Ay" 


Ar  , A0 


c 

m 


e,  e 


Temperature 
Total  temperature 

Cartesian  velocity  component,  Eq.  (1-22) 

Reference  speed,  meters/sec 
Velocity  vector 

Contravariant  velocity  component 
Physical  velocity  component 
Computat ional  curvilinear  coordinate 
Spatial  differencing  parameter,  Eq.  (1-20) 

Temporal  differencing  parameters,  Eqs.  ( I — 1 1 , 1-12) 
Specific  heat  rates,  cp/cy 
Kronecker  delta 
Time  increment 

Mesh  spacing  for  cartesian  coordinate  xm,  Eq.  (1-22) 

Mesh  spacing  for  coordinate  y® 

Mesh  spacing  for  computational  polar  coordinates  r and  6 
Spatial  difference  operators,  Eq.  (1-20) 

Convergence  parameter,  Eq.  (1-34) 

Artificial  viscosity,  Eq.  (1-22) 

Angular  polar  computational  and  physical  coordinates 
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LIST  OF  SYMBOLS  - PART  I (cont'd) 


p Dynamic  viscosity 

p Density 

4 Dependent  variable  column  vector,  Eq.  (1-9) 

4 Dependent  variable  column  vector,  Eq . (I-15a) 

a)  Fourth  derivation  dissipation  coefficient,  Eq.  (1-23) 

Subscript  s 

i,j,k  Denote  cqvariant  tensor  components  or  function  evaluated  at  grid 
point  (yi,  yj,  yk) 

r Denotes  dimensional  reference  value 

t*0  Denotes  inviscid  initial  solution 

,i  Denotes  covariant  derivative  with  respect  to  coordinate  y1 

Superscripts 

i,j,k  Denote  contr avariant  tensor  component  or  grid  point  location 

L Denotes  lower  periodic  boundary  of  blade  passage 

U Denotes  upper  periodic  boundary  of  blade  passage 

*,  **  Denote  intermediate  solutions  of  alternating  direction  procedure 
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BACKGROUND 


The  UTRC  Navier-Stokes  computer  code  obtains  implicit  time-marching 
finite  difference  solutions  to  a general  coordinate  representation  of  the 
equations  of  motion.  The  numerical  solution  procedure  on  which  this  code  is 
based  was  originally  developed  by  Briley  and  McDonald  (Ref.  1-1).  A typical 
time  step  in  their  procedure  consists  of  a time-wise  linearization  followed  by 
a fully  implicit  difference  approximation  which  is  solved  by  an  ADI  (Alternat- 
ing Direction  Implicit)  procedure  of  the  Douglas-Gunn  type  (Ref.  1-2).  The 
advantage  with  ADI  methods  is  that  a short  sequence  of  simple  matrix  inver- 
sions replaces  the  complicated  matrix  inversion  problem  associated  with  a 
direct  solution  of  the  implicit  equations.  In  this  way  a real  savings  in 
computer  time  is  made  without  sacrificing  accuracy  or  stability.  Briley, 
McDonald,  and  Gibeling  (Ref.  1-3)  have  shown  that  the  ADI  scheme  ha6  run 
stably  and  accurately  with  time  steps  which  are  orders  of  magnitude  larger 
than  the  explicit  stability  limit  (Ref.  1-4).  The  first  applications  of  the 
implicit  time-marching  Navier-Stokes  computer  code  involved  the  calculation 
of  laminar  and  turbulent  flows  in  rectangular  ducts  (Refs.  1-1  and  1-3). 
Recently  the  basic  numerical  analysis  (and  computer  code)  ha6  been  extended 
by  Gibeling,  Shamroth  and  Eiseman  to  consider  a general  curvilinear  coordi- 
nate form  of  the  Navier-Stokes  equations.  In  addition,  they  applied  the 
curvilinear  coordinate  version  of  the  Navier-Stokes  code  to  calculate  low  to 
moderate  Reynolds  number  external  flows  past  simple  aerodynamic  shapes  (Ref. 
1-5). 


In  the  present  effort  sample  external  flow  calculations  have  been  per- 
formed to  assess  the  reliability  and  accuracy  of  the  current  version  of  the 
implicit  time-marching  Navier-Stokes  computer  code,  and  this  code  ha6  been 
modified  to  treat  cascade  blade  passage  flows.  The  governing  equations  and 
numerical  solution  procedure  upon  which  the  UTRC  Navier-Stokes  computer  code 
is  based  will  be  outlined  below  for  the  special  case  of  laminar,  steady  adia- 
batic flow.  The  equations  of  motion  and  the  numerical  solution  procedure  are 
described  for  three-dimensional  flows;  however,  only  two-dimensional  applica- 
tions have  been  considered  herein.  For  a more  general  treatment  of  the  gover- 
ning equations  and  more  complete  details  on  the  numerical  method  the  reader 
i6  referred  to  References  I— 1 , 1-3,  and  1-5. 
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EQUATIONS  OF  MOTION 


All  quantities  in  the  following  equations  are  dimensionless;  physical 
velocities  are  normalized  by  Uf,  density  by  Pf , pressure  by  PrU  r,  tem- 
perature by  T , molecular  viscosity  by  Pr,  length  by  Lr,  and  time  by  Ur/Lr. 

The  fluid  is  assumed  to  be  a perfect  gas  with  zero  bulk  viscosity  coefficient 
and  constant  molecular  viscosity  and  specific  heat.  The  flow  i6  assumed  to 
be  adiabatic  with  constant  total  temperature,  T , and  body  forces  are  assumed 
to  be  negligible.  Under  these  assumptions,  the  Navier-St okes  equations  can 

1 o o 

be  written  for  time-independent,  curvilinear  coordinates  (y1 , y , yJ)  as  follows 
(Ref.  1-6):  the  continuity  equation  is 

+ (/> v* ),  j = o (i-D 


and  the  momentum  equations  are  given  by 

' (If  * v' vi' l)  • - 9"  “•)  * 55?  g"  Al>i + 


(1-2) 


where  p is  the  density,  v is  the  velocity,  p is  the  pressure,  p is  the 
molecular  viscosity,  Re  * prUrLr^ur  Reyno^8  number.  In  addition,  g1^ 

are  components  of  the  inverse  metric,  v1  are  contravariant  velocity  components, 
subscripts  after  commas  denote  covariant  derivatives,  the  indices  i,  j,  and  k 
vary  from  1 to  3,  and  repeated  indices  are  to  be  summed.  The  energy  equation 
can  be  replaced  with  an  adiabatic  equation  of  state 

p = p |a  + Bg(. v'vij  (1-3) 


where  A * 1 J yMf^ , B ■ (y  - l)/2y,  g — are  components  of  the  metric  tensor, 
Mf  is  the  reference  Mach  number,  and  y is  the  specific  heat  ratio  (Cp/cv). 

After  substituting  Eq.  (1-3)  into  Eq.  (1-2)  to  eliminate  the  pressure, 
the  governing  equations  can  be  expressed  in  a convenient  form  for  numerical 
calculations  (Ref.  1-3).  It  follows  after  some  algebra  that 


* A_  Ipjvh  *o 
41  V 


(1-4) 
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and 


♦ L 


k dv 


%km  a2vJ 


ay 


k.  m + V + 'ijV 


ay  dy 


♦ Gijk  ',v‘vk  =° 


(1-5) 


The  coefficients  in  Eq . (1-5)  are  defined  by 
K = AJ 

p rn  / a^rn  \ 

F|k,  !K!l  ♦•9|k*i)  J 

G . - J.  i252.  /Sm."  „„  _mn\ 

i,m  2 ay'  n VeV9  ) 


R.  = A 


At  «J  _ J . Jk 


ay  • - ay 

Dum  * w j [-§-»r»r-««^"-*r»T] 

aSmk 


, = iiil  + _S_ 

'»  dyk  2Re 


. mk 


ay' 


. nm  pk  ^9 

jg  gK 


EH  _L  _£J_  nmk 

ayj  J ay'  9 . 


(1-6) 


, k k ^Gij_  fi  ^Qnm  f„kn*rn  km,n  2 nm ftk"| 

L'i'E'i  + l/F-t2W-’l^[5  V Vl9  ®ij 

where  the  Kronecker  symbol  {^vanishes  unless  i * j in  which  case  it  is 
unity  and 


k - JL  [-S-  *k  _ mk  ^9jm 

ij  - »«  [7  4yi  »,  ag  4>J 


Equations  (1-4)  and  (1-5)  are  the  Navier-Stokes  equations  in  a fixed  coordinate 
frame  (y1 , y2 , y^)  with  the  density,  p,  and  the  contravariant  velocity  compo- 
nents, v1,  i • 1,  2,  3,  taken  as  dependent  variables.  Once  the  curvilinear 
coordinate  system  is  prescribed,  the  Jacobian  and  metric  tensor  components 
and  hence  the  coefficients,  defined  by  Eqs.  (1-6)  and  (1-7),  become  know 
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functions  of  the  spatial  coordinates.  Previous  experience  (Ref.  1-5)  has 
indicated  that  the  use  of  contravariant  velocity  components  as  dependent 
variables  can  lead  to  serious  truncation  errors  and  as  a result  the  current 
version  of  the  UTRC  Navier-Stokes  computer  code  has  been  formulated  to 
consider  the  density  and  the  physical  velocity  components,  v(i)  as  basic 
dependent  variables,  where 


v(i)  = (9ji) 


1/2 


(no  sum  on  i)  (1-8) 


The  appropriate  form  of  the  Navier-Stokes  equations  is  then  obtained  after 
substituting  Eq . (1-8)  into  Eqs.  (1-4)  and  (1-5). 

Although  the  use  of  generalized  (nonorthogonal ) coordinates  leads  to  a 
complicated  form  of  the  equations  of  motion.  Such  coordinates  offer  signifi- 
cant computational  advantages.  In  particular,  physical  boundaries  of  the 
flow  region  can  be  represented  as  coordinate  surfaces,  thus  removing  the 
need  for  fractional  cells  and  boundary  condition  interpolations.  Further, 
a uniform  mesh  can  be  used  in  computational  space  and  mapped  into  a mesh 
suitably  distributed  in  physical  space  to  capture  large  solution  gradients 
such  as  those  occurring  in  boudary  layers  and  near  airfoil  leading  and 
trailing  edges.  Finally,  the  uniform  mesh  in  computational  space  simplifies 
the  finite  difference  approximations  of  derivative  terms. 
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NUMERICAL  SOLUTION  PROCEDURE 


The  numerical  method  can  be  briefly  outlined  as  follows:  the  governing 
equations  are  replaced  by  an  implicit  time  difference  approximation,  option- 
ally a backward  difference  or  Crank-Nicolson  scheme.  Terms  involving  non- 
linearities  at  the  implicit  time  level  are  linearized  by  Taylor  expansion  about 
the  solution  at  the  known  time  level,  and  spatial  difference  approximations  are 
introduced.  The  result  is  a system  of  linear  difference  equations  for  the 
dependent  variables  at  the  unknown  or  implicit  time  level.  To  solve  these  dif- 
ference equations,  the  Douglas-Gunn  (Ref.  1-2)  procedure  for  generating 
alternating-direction  implicit  (ADI)  schemes  is  introduced.  This  technique 
leads  to  systems  of  coupled  linear  difference  equations  having  narrow  block- 
banded  matrix  structures  which  can  be  solved  efficiently  by  standard  block- 
elimination  methods. 

To  describe  the  numerical  procedure  it  is  convenient  to  express  the 
Navier-Stokes  equations  in  the  following  matrix  form: 

W + s(<#>)  d-9) 

0 1 


where  $ is  a column  vector  containing  the  dependent  variables,  H is  a column 
vector  function  of  if,  is  a column  vector  whose  elements  are  functions  of  the 
dependent  variables  and  their  spatial  derivatives  in  a single  coordinate  direc- 
tion, and  S is  a column  vector  whose  elements  are  functions  of  the  mixed 
second  order  spatial  derivatives  of  the  dependent  variables.  For  example, 
when  the  contravar iant  velocity  components  are  treated  as  dependent  variables, 
$,  H ($)  and  S ($)  have  the  form  (cf.  Eqs.  (1-4)  and  (1-5)): 


\ip 

p 

V1 

V2 

H = 

JQij  P*] 

V3 

,Jfl3  ,'v'. 

S = 


av 

ay'dyS 


MD 


V «) 


K * 0 S + K * ■>?;) 


(1-10) 
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and  the  elements  of  the  column  vector  ^(if)  consist  of  the  remaining  terms  in 
Eqs . (1-4)  and  (1-5).  It  should  be  noted,  that  certain  definitions  and  descrip- 
tions of  terms  provided  in  the  present  discussion  differ  from  those  given 
earlier  in  Refs.  1-1,  1-3  and  1-5.  However,  the  definitions  used  here  are 
appropriate  to  the  current  version  of  the  implicit  time-marching  Navier-Stokes 
computer  code. 


The  solution  domain  is  discretized  by  grid  points  having  equal  spacings 
in  the  computational  coordinates,  Ayl,  Ay^ , and  Ay  in  the  y^,  y and  y direc- 
tions, respectively,  and  an  arbitrary  time  step.  At.  The  subscripts  i,  j,  k 
and  superscript  n are  grid  point  indices  associated  with  y^,  y^,  y^  and  t, 
respectively,  and  thus  ^ denotes  $ (y^,  y?,  y^,  tn) . It  is  assumed  that 

the  solution  is  known  at  ink  n level,  tn,  and  is  desired  at  the  (n  + 1)  level, 
tn+1.  At  the  risk  f an  occasional  ambiguity,  one  or  more  of  the  subscripts 
is  frequently  omitted,  so  that  4>n  is  equivalent  to  j 


Linearization  Scheme 

A linear  difference  approximation  to  the  nonlinear  governing  equations  is 
obtained  from  the  following  time-difference  replacement  of  Eq.  (1-9): 

(H"4,-Hn)/At  =P,  ^"*'+0 -0,)2>n  + /3*  Sn  + ,+  (l-/32)Sn  (i-ll) 

where,  for  example,  Hn4^  « H(4>n+^).  The  parameters  (0  £ 6^,  8 £ 1)  permit 

variable  centering  of  the  scheme  in  time.  Thus  Eq.  (I-ll)  produces  a backward 
difference  or  fully  implicit  formulation  for  6^  * 62  “ lt  a Crank-Nicolson 
formulation  for  8^  - B2  * 1/2,  and  a forward  difference  or  fully  explicit 
scheme  for  * 62  * Unconditional  stability  is  anticipated  for 
®1*  > 1/2*  In  the  present  method  the  column  vector  S ($)  is  treated 

explicitly.  Thus  with  8 ■ 8^,  Eq.  (I-ll)  reduces  to 

(Hn4,-Hn)/At«0,  V + ' + (l-/5)  V+S°  (1-12) 

The  mixed  derivative  terms  (in  S($))  could  be  treated  implicitly  within  the 
ADI  framework;  however,  this  would  increase  the  number  of  intermediate  steps 
and  thereby  complicate  the  solution  procedure.  Test  cases  computed  while 
developing  the  present  numerical  method  (Refs.  1-1  and  1-3)  have  indicated 
that  the  explicit  treatment  of  mixed  derivative  terms  has  no  observable  adverse 
affect  on  stability. 
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The  linearization  is  accomplished  by  a two  step  process  of  expansion 
about  the  known  time  level  tn  and  subsequent  approximation  of  the  quantity 
(3<f>/3t)n  At,  which  arises  from  chain  rule  differentiation,  by  (4>n+-*  - $n) , 

The  result  is 

Hn*'s  H°+  (dH/d*)n(*n+'-*n)  +0  (At)2  (I-13a) 


H=  3%-  (d%/d4>f 


( I-13b) 


The  matrix  3H/34>  is  a standard  Jacobian  whose  elements  are  defined  by 
(3H/34))qr  = 3Hq/3if>r.  The  operator  elements  of  the  matrix  3.2>/34>  are  similarly 
ordered,  i.e.,  (3^/3if)qr  « 3^q/3<t>r  however,  the  intended  meaning  of  the 
operator  elements  requires  some  clarification.  For  the  qth  row.  the  operation 
(3  2kq/3<t>)n  (<f>n+^  - 4>n)  is  understood  to  mean  that  { 3 / 3 1 2) > [<Ky  , y^,  y , t)]  }nAt 

is  computed  and  that  all  occurrences  of  (3<f>r/3t)n  arising  from  chain  rule 
differentiation  are  replaced  by  (4>rn+l  - <J>£  )/At.  The  substitution  of  Eqs. 

(1-13)  into  Eq.  (1-12)  leads  to  the  following  linear,  implicit,  first  order 
accurate  time-differenced  scheme: 


+ ♦ sn 


(1-114) 


Equation  (1-14)  is  linear  in  the  quantity  (4>n+^  - 4>n)  and  all  other  quantities 
are  either  known  or  evaluated  at  the  nth  time  level.  It  is  convenient  to 
solve  Eq.  (1-14)  for  (<t>n+^  - $n)  rather  than  (fin+1.  This  reduces  roundoff 
errors,  since  it  is  presumably  better  to  compute  a small  0 (At)  change  in  an 
0 (1)  quantity  than  the  quantity  itself.  After  defining  the  symbols: 


* ■ 

(I-15a) 

-for 

(I-15b) 

(dS/d+f 

(I-15c) 

Equations  (1-14)  can  be  written  in  the  following  simplified  form 

(*  + £k\j)  *n  + '«  At[3>n+Sn]  (1-16) 


1 
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Alternating  Direction  Procedure 


The  solution  of  Eq.  (1-16)  is  accomplished  by  application  of  an  alter- 
nating-direction implicit  (ADI)  technique  which  is  a generalization  of  the 
procedure  developed  by  Douglas  and  Gunn  (Ref.  1-2)  for  generating  ADI  schemes 
as  perturbations  of  fundamental  implicit  difference  schemes  such  as  the  back- 
ward-differences or  Crank-Nicolson  schemes.  The  vector  operator  2(<P)  contains 
terms  which  are  fund  ions  of  <p  and,  in  addition,  terms  which  are  functions  of 
if  and  the  first  and  second  order  derivatives  of  if  with  respect  to  y^,  , and 

y^,  but  no  mixed  derivatives.  Thus  2 (and  hence  J)  can  be  split  into  three 
operators,  associated  with  the  y^-,  y^,  and  y^  coordinates,  each 

having  the  functional  formJ^m  = ^(<f,  3/8ym,  a^/Sy^Sy1”)  . Those  terms  of 
J^(4>)  which  do  not  contain  a spatial  derivative  are  grouped  under  the  operator 
Equation  (1-16)  then  becomes 

[A  + At(^|  + ^2+-^3)]^n  + l=At  [jfc"  +Snj  (1-17) 

and  the  Douglas-Gunn  representation  of  Eq . (1-17)  can  be  written  as  the  fol- 
lowing three  step  solution  procedure: 

(A  + At/)V'*  = At[X  + ^2  + ^3  + 5n]  (I-18a) 

(A  ♦ (I-18b) 

(A  -6  At-/3)  ^n  + l«  A***  (I-18c) 


where  if*  and  if**  are  intermediate  solutions.  If  if*  and  if**  are  eliminated, 
Eqs.  (1-18)  become 

(a  + At/)AH  (a  + At-/2)AH  (A  + At^3)tn+'=  At  + 3>2  + 2>*  ♦ Sn]  (1-19) 

and  after  performing  the  multiplication  on  the  left-hand  side  of  Eq.  (1-19), 
it  is  apparent  that  Eq.  (1-19)  approximates  Eq . (1-17)  to  order  (A t)^. 

Spatial  differencing  of  Eqs.  (1-18)  is  accomplished  by  replacing  deri- 
vative operators  such  as  9/8yID  and  8^/3ym9yln  (no  sum  on  m)  by  corresponding 
three  point  finite  difference  operators,  D^  and  D^  where 
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( I-20a) 


D m4>=  [aA_+ (l-a)A+j^/Ayrn  = d<t>/dyW  + (?  [(Ay)2  + (a  - l/2)Aym] 

4>  = (A+A_)  </>/(Aym)2=  dZ4>/dyT'  dy™  + 0 (Aym)  (I-20b) 


Here,  £+4>  and  A_d>  represent  forward  and  backward  differences  in  (}; . For 
example,  the  difference  approximation  to  the  derivative  of  4>  at  the  point 
(y}>  y j , y^)  in  the  y^  - direction  (m  = 2)  is  obtained  by  setting 

A+*i,j,k  - ♦i.j+l.k  ‘ *i,J,k  and  A_  ^ i , j , k = * i, j ,k  " ♦ i , j -1 ,k-  The 
parameter  a has  been  introduced  (0  <_  a 1)  in  Eq.  (1-20)  to  permit  con- 
tinuous variation  from  backward  to  forward  differences.  The  standard  central 
difference  formula  is  recovered  for  a = 1/2  and  was  used  for  the  numerical 

calculations  reported  here. 

2 

With  the  introduction  of  the  spatial  difference  operators  D^  and  D~, 
defined  in  Eqs.  (1-20),  the  solution  procedure  for  the  alternating  direction, 
implicit  time  marching  form  of  the  governing  equations,  Eqs.  (1-18),  can 
be  described.  Since  Dm  and  are  three  point  difference  operators,  the 
finite  difference  approximation  to  Eq.  (I-18a)  contains  i_ljj  ,k»  ^*i  j k’ 
and  j k as  unknowns.  Hence,  the  system  of  linear  equations  generated 

by  writing  Eq.  (I-18a)  at  successive  grid  points  (y^,  y ? , y^;  i = 1,  . . . , I) 
can  be  written  in  block-tridiagonal  form;  i.e., 


+ b 


n 


* 


* 

i + l,j,k 


(1-21) 


where  a,  b,  and  c are  square  matrices  and  d is  a column  vector,  each  con- 
taining only  n-level  quantities.  When  applied  at  successive  grid  points 

(i  = 1,  , I ),  Eq.  (21)  generates  a block-tridiagonal  system  of  equations 

for  ip*  which,  after  appropriate  treatment  of  boundary  conditions,  can  be 
solved  efficiently  by  using  standard  block-elimination  methods  (Ref.  1-7). 

The  solution  procedure  for  Eqs.  (I-18b,  c)  is  analogous  to  that  just  described 
for  Eq.  (I-18a)  . 


Artificial  Dissipation 

In  computing  solutions  for  high  Reynolds  number  flows,  it  is  often 
necessary  to  add  a form  of  artificial  viscosity  or  dissipation.  One  possible 
dissipation  term  in  common  use  is  based  on  an  observation  by  Roache  (Ref.  1-8) 
that  for  a linear  model  problem  representing  a one-dimensional  balance  of 
convection  and  diffusion  terms,  solutions  obtained  using  central  differences 
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for  the  convection  terms  are  well  behaved  provided  the  mesh  Reynolds  number 
Re^xtn  * jum!Axm  Re  is  2,  but  that  qualitative  inaccuracies  associated 
with  boundary  conditions  occur  for  Re^xm  > 2.  This  suggests  the  use  of  an 
artificial  viscosity  term  of  the  form  emD2$ , where 


_ 


,|a> 


i 

Re 


i 

Re 


04 


ReAv  > 2 


°'ReA«m-2 


(1-22) 


to  insure  that  the  local  effective  mesh  Reynolds  number  is  no  greater  than  two. 
This  result  has  been  extended  for  generalized  tensor  equations  and  the  appro- 
priate dissipation  terms  for  the  continuity  and  momentum  equations  have  been 
incorporated  by  Gibeling  et  al.  (Ref.  1-5)  into  the  present  version  of  the 
UTRC  Navier-Stokes  computer  code. 


A second  type  of  artificial  damping  which  is  a fourth-order  dissipation 
term  has  been  suggested  by  Beam  and  Warming  (Ref.  1-9)  to  damp  small  wave- 
length disturbances.  In  the  present  formulation  an  explicit  fourth-order 
damping  term  was  added  directly  to  the  fundamental  difference  scheme, 

Eq.  (1-16),  as  follows: 


(A  + At  J)  + ' = At  [.2>  l$n) 


1 8 


aV 

d(ymf 


(1-23) 


This  dissipation 
matrix  structure 


term  is  treated  explicitly  to  retain  the  block  tridiagonal 
of  the  finite  difference  form  of  the  governing  equations. 
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NUMERICAL  RESULTS  FOR  FLOW  PAST  A CIRCULAR  CYLINDER 


Under  the  present  effort  implicit  time-marching  Navier-Stokes  solutions 
have  been  computed  for  symmetric,  two-dimensional  flows,  at  Reynolds  numbers 
(based  on  diameter)  of  forty  and  eighty,  past  a circular  cylinder  (Fig.  1). 

In  both  cases  the  free  stream  Mach  number  was  set  equal  to  0.2.  The  circular 
cylinder  case  is  a convenient  example  because  it  represents  a relatively 
simple  geometry  and  both  experimental  data  and  other  numerical  solutions  are 
available  for  comparison  (cf.  Ref.  1-10) . In  particular,  Gibeling,  Shamroth, 
and  Eiseman  (Ref.  1-5)  have  computed  the  flow  past  a circular  cylinder  at  a 
Reynolds  number  of  forty  using  the  present  numerical  method.  These  authors 
considered  the  density  and  contravar iant  velocity  components  as  dependent 
variables  and  did  not  use  artificial  fourth-order  dissipation  (cf.  Eq.  (1-23)) 
in  their  calculations.  However,  they  later  modified  the  Navier-Stokes  computer 
code  to  consider  density  and  the  physical  velocity  components  as  dependent 
variables  and  introduced  fourth  order  dissipation  into  the  calculation  proce- 
dure in  order  to  compute  flows  past  a symmetric  Joukowski  airfoil  (Ref.  1-5). 
Circular  cylinder  calculations  are  repeated  here  both  to  test  the  current 
version  of  the  Navier-Stokes  computer  code,  and  so  that  the  present  investi- 
gators could  gain  experience  in  working  with  the  code  by  applying  it  to  a 
simple  external  flow  configuration. 

Due  to  symmetry  results  are  computed  for  only  the  upper  half-plane  of 
the  flow  field  using  a 35  x 35  mesh  embedded  in  a polar  coordinate  system. 

The  outer  boundary  of  the  computational  region  was  taken  to  be  fifteen 
diameters  from  the  cylinder  center.  A nonuniform  grid  spacing  in  physical 
space  is  applied  with  points  along  radial  lines  concentrated  near  the  cylinder 
boundary  and  points  on  the  circumferential  lines  concentrated  near  the  front 
and  rear  stagnation  points.  Mesh  points  are  distributed  according  to  the 
relations 


8(iA0)  = 7r/2 -ir/2  tanh  { Dg  [ I -2  ( i-l)  A0/tt  j }/tanh  Dg  , i = I .....  ,1  (I-24a) 


R(  jAr)=  30-29  tanh  {Dr  [l-(  j-l)Ar]  } / tanh  Dr 


i * i .....  j 


( I — 2 A b) 


where  A9  * ir/(I-l),  Ar  = 1/(J-1),  Dg  = 0.75,  and  Dr  = 2.7  are 
damping  constants  for  the  0 and  R directions,  respectively,  and  I « J * 35. 
The  coordinate  distributions,  Eqs.  (1-24),  are  special  cases  of  the  Roberts 
boundary  layer  transformation  (Ref.  1-11)  and  are  shown  plotted  in  Fig.  2. 
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The  potential  flow  past  a circular  cylinder  with  a compressibility 
correction  on  density  is  specified  as  the  initial  flow  field  for  the  present 
viscous  time  marching  solutions  with  the  exception  that  the  initial  velocity 
on  the  surface  of  the  cylinder  is  set  equal  to  zero.  Thus,  the  initial  condi- 
tions are  as  follows: 


( '-rMr2/4)/2  j/[MMr  v)2  (l-y(Mr  v)2  /4)/z] 

v(i ) | =-( i«  R'2 ) sin  8,  R*l 

= 0,  R = I 
v (2)|  = ( I * R 2 ) cos  e 


( I-25a) 

( 1—2  5b) 

(I-25c) 


where  v(l)  and  v(2)  are  the  physical  velocity  components  in  the  R and  0 - 
directions,  respectively  and  v is  the  magnitude  of  the  velocity.  The  implicit 
time  marching  solutions  are  constrained  to  satisfy  the  following  steady 
boundary  conditions.  Symmetry  conditions  are  applied  on  the  stagnation 
streamlines;  i.e., 


dp/d&  = \/(\)--d  v(2)/d0:O,  on  0 =0,  ir  (1-26) 


thus  permitting  calculations  to  be  restricted  to  the  upper  half-plane  of  the 
flow.  The  velocity  and  normal  pressure  gradient,  3p/3n,  is  set  equal  to  zero 
at  the  cylinder  surface.  This  condition  is  required  to  determine  a value  of 
density  at  the  surface  and  in  the  present  case  of  constant  total  temperature, 
this  condition  is  equivalent  to  3p/3n  = 0.  Therefore 


dp/d R = v(l)=  v(2)-0,  on  R = l 


(1-27) 


Finally,  velocity  and  density  are  set  to  their  inviscid  values  over  the 
upstream  three  quarters  of  the  outer  boundary  while  first  derivatives  of  the 
physical  velocity  components  are  set  to  zero  and  the  pressure  is  set  to  its 
inviscid  value  over  the  remainder  of  the  far-field  boundary;  i.e., 

I t = 0 ,V(l)t  V(l>lt:0  •V(2)lV(2)lt*0  ' (I’28) 


on  R = 30,  tt / 4 < 0<  w 
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and 


d v ( I ) / d R = <Jv(2)  / dR  --0, 


p--p  t _ / {a-*-b[v2(i)  + v2(2)]|,  on  R = 30,0<e<W4 


(1-29) 


The  last  condition  follows  from  Eqs.  (1-3)  and  (1-8)  with  = 0 due  to 
orthogonality. 

At  this  point  it  is  convenient  to  review  the  solution  procedure  using  the 
circular  cylinder  problem  as  a vehicle  to  illustrate  the  method  of  incorporat- 
ing the  boundary  conditions  into  the  finite  difference  approximation.  The 
foregoing  steady  state  boundary  conditions  on  the  lines  G =0  and  it  and 
R = 1 and  30  can  each  be  expressed  in  the  form 

(1-30) 


where  is  a matrix  operator  containing  derivatives  in  the  coordinate 

directions  and  Sn  is  a column  vector  function  of  <t>n.  Hence  the  intermediate 
solution,  after  the  nth  time  step  of  the  alternating  direction  implicit 

procedure,  is  calculated  on  the  circumferential  lines,  R (jAr)  = constant, 
j = 2,  ...,  J-l,  as  a solution  of  the  equation  set 


(I-31a) 

(A*  At  /,)|tj  **  J +fi2n*S"]itj  . i =2... 

...r-1 

(l-31b) 

B 1 .0,*  =J3  !,.<#>" 

BC  1 t,j  *|,j  °BC  1 1.1  ^ 

(l-31c) 

After  approximating  derivatives  at  boundaries  by  three  point  finite  difference 
expressions,  a block  tridiagonal  system  of  linear  equations  of  the  form 
(cf.  Eq.  (1-21)) 
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. n ,*  n , * n . * .n 

bl,  j ti,  i ♦c  l.l  ^2,  i + Q l.J  ^3,  i = d* , j 


( I-32a) 


n .*  . n ,*  n .*  n 

°i,j  '/'i-l.j*  b,,j  't'i.i  4Ci,j  --djj  .1=2 I -I  ( I — 32b) 


Cn  a,*  + a",  u,*  ♦ b":  **  sd", 

l.i  M-2.j  J.)  M-I.i  M.J  1,1 


(I-32c) 


results.  Values  of  are  determined  by  using  Gaussian  elimination  methods  to 
solve  this  system  of  equations.  Similarly,  the  viscous  flow  field  at  the 
(n  + 1)  th  time  level,  , is  obtained  as  a solution  of  the  equations 


^BC  lj,l 

d,"4' 

r i.i 

= -fiec 

<*>"  ♦ Sn 

I i,i  ^ i.i  i.i 

(I-33a) 

(A.  At  *2)u 

n*i 

. * 

* hi 

1 A',i 

Ki  . j = 2. ..  . j- i 

(I-33b) 

*BC  li.J 

yL 

y i,J 

= Bbc 

^ i.j 

(I-33c) 

on  the  radial  lines,  G(iA0)  = 

constant,  i 

*=2,  ...  1-1,  after 

spatial  dif- 

ference  approximations  are  introduced  to  express  the  foregoing  equations  in 
block  tridiagonal  form.  The  viscous  solution  evolves  from  the  inviscid  flow 
field,  i|°,  until  the  difference  between  the  flows  at  two  successive  time  steps, 
N and  N + 1,  satisfies  an  imposed  convergence  criterion.  Specifically,  the 
maximum  value  over  the  entire  field  of  the  absolute  difference  between  the 
dependent  variables  at  successive  time  steps  must  be  less  than  some  small 
number , c ; i. e. , 


N + i N 
<p  - <p  \\  -. 


nh 


N 


-*k 


(1-34) 


Here  the  subscript  k refers  to  the  kth  component  of  the  column  vector  $ • In 
the  solution  procedure  the  magnitude  of  the  time  step  at  the  nth  time  level, 
n » 1 , . . . , N , is  a variable  parameter  which  depends  on  the  value  of 
I I _ *n-l | I 


In  the  present  study  implicit  time-marching  Navier-Stokes  solutions 
were  advanced  through  160  and  200  time  steps  for  circular  cylinder  flows  at 
Reynolds  numbers  of  40  and  80,  respectively.  The  dimensionless  time  step  used 
in  these  calculations  varied  between  0.15  and  0.5.  Artificial  fourth  order 
dissipation  was  at  first  applied  for  the  circular  cylinder  calculations,  but 
was  found  to  lead  to  diverging  solutions  and  spurious  results  In  the  far  field. 
Hence,  the  only  artificial  viscosity  term  used  below  was  that  based  on  the  mesh 
Reynolds  number  criterion,  Eq.  (1-22).  The  computations  were  carried  out  on 
the  UTRC  Univac  1110  computer  system  and  the  computing  time  for  the  nonorthog- 
onal  form  of  the  governing  equations  used  for  the  circular  cylinder  flows  is 
approximately  0.95  CPU  minutes  per  time  step  or  7.8  x 10”^  CPU  minutes  per 
grid  point  per  time  step.  Selected  results  from  the  present  calculations  for 
flow  past  a circular  cylinder  are  described  below. 

The  variation  of  the  minimum  pressure  and  the  separation  angle,  0fi,  with 
time  is  depicted  in  Fig.  3 for  the  Re  = 40  case  and  in  Fig.  4 for  the 
Re  = 80  case.  After  160  time  steps  the  separation  angle  reaches  a value  of 
51.4°  for  Re  = 40.  This  is  to  be  compared  with  values  of  50.0°,  52.5°,  53.7°, 
and  53.9°  obtained  by  Kawaguti  (Ref.  1-12),  Apelt  (Ref.  1-13),  Kawagutl  and 
Jain  (Ref.  1-14)  and  Son  and  Hanratty  (Ref.  1-15).  The  development  of  the 
pressure  distribution  around  the  surface  of  the  cylinder  with  time  is  shown 
in  Fig.  5 for  Re  = 40  and  in  Fig.  6 for  Re  = 80.  The  time  history  plots. 

Figs.  3 through  6 reveal  that  although  the  implicit  time-marching  results 
appear  to  be  converging,  steady  state  solutions  have  not  been  achieved  even 
after  a considerable  number  of  iterations  or  time  steps.  The  maximum  dif- 
ference between  dependent  variables  at  successive  time  steps  (cf.  Eq.  (34)) 
is  of  the  order  of  2 x 10-^  after  160  (Re  = 40)  and  200  (Re  - 80)  timA  steps. 

In  addition,  as  will  be  seen  below,  the  results  after  160  time  steps  for  the 
Reynolds  number  of  forty  case  are  not  in  good  agreement  with  previous  cal- 
culations or  experimental  measurement. 

The  present  prediction  of  surface  pressure  distribution  for  Re  * 40 
after  160  time  steps  is  compared  with  the  predictions  of  Son  and  Hanratty 
(Ref.  1-15),  Kawaguti  (Ref.  1-12)  and  Gibeling,  Shamroth,  and  Eiseman 
(Ref.  1-5)  in  Fig.  7.  Note  that  the  angular  coordinate,  0,  is  equal  to  180° 
at  the  front  stagnation  point.  The  Kawaguti  pressure  coefficient  prediction 
was  given  relative  to  the  rear  stagnation  point  pressure  and  in  Fig.  7 the 
Kawaguti  rear  stagnation  point  pressure  was  arbitrarily  set  at  the  Son  and 
Hanratty  value.  The  predictions  of  Refs.  (1-15)  and  (1-12)  were  both  obtained 
from  solutions  of  the  incompressible  Navier-Stokes  equations  while  the  solution 
of  Ref.  1-5  and  the  present  solution  were  obtained  from  the  compressible 
equations  with  Mach  number  equal  to  0.2.  As  is  clear  from  Fig.  7,  the  agree- 
ment between  the  present  predictions  and  those  of  Refs.  1-5,  1-12,  and  1-15 
is  rather  poor  with  the  present  analysis  over  predicting  the  pressure  over 
most  of  cylinder  surface.  In  addition,  the  time  history  plots,  Figs.  3 and 
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5,  indicate  that  with  continued  iteration  the  pressure  prediction  based  on  the 
present  procedure  would  deviate  still  further  from  previous  results.  Indeed, 
the  results  achieved  with  the  present  method  after  80  time  steps  are  in  better 
agreement  with  the  previous  predictions  than  those  achieved  after  160  time 
steps.  Wake  centerline  velocity  distributions  for  Reynolds  numbers  of  40  and 
80,  based  on  the  present  analysis,  along  with  the  predictions  of  Kawaguti 
(Ref.  1-12)  and  Gibeling  et  al  (Ref.  1-5),  and  the  experimental  measurements 
of  Coutanceau  and  Bouard  (Ref.  1-10)  for  a Reynolds  number  of  40  are  depicted 
in  Fig.  8.  Again  the  agreement  between  the  data  and  the  present  results  for 
Re  * 40  is  not  good.  As  a final  illustration  of  the  results  of  the  present 
calculations,  the  velocity  profiles  for  Re  * 40  and  Re  = 80  at  several  azimuthal 
locations  are  shown  in  Fig.  9 along  with  the  predictions  of  Gibeling,  Shamroth, 
and  Eiseman  (Ref.  1-5)  for  Re  = 40.  Substantial  differences  can  be  observed 
between  the  results  of  these  two  analyses. 

Based  on  the  present  study  the  performance  of  the  implicit  time-marching 
Navier-Stokes  computer  code  has  not  been  satisfactory.  Converged  solutions 
could  not  be  obtained  for  simple  low  Reynolds  number  flows  and  the  results 
after  a considerable  number  of  time  steps  are  not  in  good  agreement  with  those 
of  previous  analyses  or  experimental  data.  In  contrast,  Gibeling,  Shamroth 
and  Eiseman  (Ref.  1-5)  used  essentially  the  same  version  of  the  Navier-Stokes 
computer  code  and  reported  a converged  solution  for  the  flow  past  a circular 
cylinder  at  Re  * 40  after  only  80  time  steps,  and  their  results  are  in  very 
good  agreement  with  the  previous  incompressible  analyses  of  Refs.  1-12  and 
1-15  and  the  experimental  data  of  Ref.  1-10.  The  major  difference  between 
the  present  analysis  and  that  of  Ref.  1-5  is  that  physical  velocity  components 
are  treated  as  the  dependent  variables  here,  while  contravariant  velocity 
components  were  assumed  as  dependent  variables  in  the  former  analysis.  The 
modifications  to  the  Navier-Stokes  computer  code  to  accomplish  this  change  of 
dependent  variables  were  made  by  Gibeling  et  al  in  order  to  decrease  numerical 
truncation  errors.  At  this  time  the  reason  for  the  discrepancies  between  the 
present  results  and  those  of  the  Ref.  1-5  study  are  not  apparent,  but  further 
studies  to  clarify  this  situation  seem  warranted.  In  this  regard  information 
concerning  the  time  step  distribution  and  convergence  criterion  used  in  the 
previous  analysis  (Ref.  1-5)  as  well  as  intermediate  results  for  pressure  and 
velocity  would  be  most  useful. 
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TWO-DIMENSIONAL  CASCADE  FLOWS 


A major  objective  of  the  present  study  has  been  to  further  develop  the 
implicit  time-marching  Navier-Stokes  computer  code  so  that  this  code  could 
be  applied  to  the  calculation  of  viscous  flow  through  cascades.  In  par- 
ticular, a solution  scheme  for  cascade  flow  fields  has  been  formulated  and 
the  relevant  subroutines  of  the  Navier-Stokes  computer  code  have  been 
modified  to  accommodate  this  scheme.  The  solution  formulation  is  based  on 
the  use  of  the  two-dimensional,  curvilinear  cascade  coordinate  system 
developed  by  Eiseman  (Ref.  1-16).  This  system  consists  of  coordinate  loops 
surrounding  each  blade  and  radial  coordinate  lines  normal  to  the  blade 
surface.  The  outermost  loop  is  constructed  so  that  cascade  periodicity 
conditions  can  be  applied  without  interpolation  between  grid  points.  The 
coordinates  are  orthogonal  on  the  airfoil  surface,  but  gradually  become 
nonorthogonal  away  from  the  airfoil.  Mesh  points  can  be  packed  in  regions 
of  large  solution  gradients  and  little  restriction  is  placed  on  airfoil  cam- 
ber and  spacing.  At  present,  Eiseman' s coordinate  generator  can  accurately 
produce  systems  for  inviscid  studies;  however,  further  refinements  are 
required  to  produce  the  higher  order  smoothness  necessary  for  viscous  cas- 
cade analyses.  An  example  of  an  Eiseman  cascade  coordinate  system  is 
depicted  in  Fig.  10. 

Consider  two  dimensional  flow  past  an  isolated  airfoil  or  a blade  in 
cascade  and  a coordinate  system  consisting  of  coordinate  loops  surrounding 
the  airfoil  (or  blade)  and  radial  coordinate  lines  emanating  from  the 
airfoil  (or  blade)  surface  and  terminating  on  an  outermost  coordinate  loop. 
The  sketch  in  Fig.  11  will  be  used  to  illustrate  the  computational  procedure 
for  both  the  flow  past  an  isolated  airfoil  and  a blade  in  cascade.  In  the 
former  case  the  reader  should  envision  the  outermost  coordinate  loop  ABCDEFGH 
in  Fig.  11  as  lying  in  the  far  field  while  in  latter,  he  should  envision  the 
outermost  loop  as  consisting  of  cascade  periodic  boundaries,  CD  and  GH,  and 
front,  HC,  and  rear,  DG,  endcaps  which  lie  in  the  far  field.  In  both  cases 
the  innermost  coordinate  loop  coincides  with  the  airfoil  (or  blade  surface). 

* 

For  the  isolated  airfoil  an  intermediate  solution,  ^ , of  the  ADI 
form  of  the  Navier-Stokes  equations,  between  the  n and  n + 1 time  levels 
is  obtained  by  solving  the  direction  - 1 equation,  i.e.,  Eq.  (I-18a)  with 
- 0 or  Eq.  (1-31),  along  the  coordinate  loops  with  the  exception  of  the 
loop  which  coincides  with  the  airfoil  surface  and  the  outermost  or  far- 
field  loop.  The  starting  radial  line  for  the  direction  - 1 calculation  can 
be  chosen  arbitrarily.  The  solution  at  the  n + 1 time  level,  ^n+^,  is  then 
obtained  by  solving  the  direction  - 2 equation;  i.e.,  Eq.  (I-I8b)  with 
4,**  . 4n+1  or  Eq.  (l-33b) , of  the  ADI  representation  on  each  radial  line 
subject  to  the  appropriate  boundary  conditions  at  the  airfoil  surface  and 
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at  the  far-field  boundary;  e.g.,  the  conditions  at  R = 1 (Eqs.  (1-27))  and 
R = 30  (Eqs.  (1-28  and  1-29))  for  the  circular  cylinder  problem.  Blade 
passage  flows  can  be  calculated  in  a similar  manner,  but  conditions  must  be 
imposed  at  the  periodic  boundaries  (CD  and  GH  in  Fig.  11)  in  lieu  of 
far-field  conditions.  These  conditions  are  that  the  density  and  velocity 
must  be  equal  at  periodic  points  (e.g.,  C and  H,  and  D and  G in  Fig.  11)  and 
values  of  these  flow  variables  on  the  periodic  boundaries  must  be  solutions 
of  the  Navier-Stokes  equations.  The  blade-to-blade  periodicit^requirement 
indicates  that  solutions  of  the  direction  - 2 equations  for  ^ should 
proceed  simultaneously  along  radial  lines  which  terminate  at  two  periodic 
points  on  the  outermost  coordinate  loop. 

If  i = 1,  2,  denotes  the  natural  basis  of  tangent  vectors  to  the 

loopwise  and  radial  coordinate  curves,  then 


(I-35a) 


=9j  j 


(I-35b) 


where  the  superscript  in  Eq.  (I-35a)  refers  to  the  contravar iant  component  of 
the  base  vector  eV.v.  The  components  of  the  unit  outward  normal  vector  to 
the  loopwise-coordinate  lines  are  obtained  by  solving  the  equations 


and 


n • 


n'  *8,2  n2  =0 


(I -3 6a) 


II  rf  11  = g j j n'  ni 


(I-36b) 


It  follows  that 


(I-37a) 


n 


2 . 


(I-37b) 


Then  since  the  velocity  must  be  equal  at  upper  and  lower  periodic  points 


27 


( V . ?(l))u  * ( v'gM  * v2  g2|  )u  = - (v1  g(|  +v2g21  )L  « -(  v-  e(1)  )L  (i-38a) 

( v • rT)^*  s ( J v 2/Vg^)U  = - ( J v2/ ■/  g , , ) = - ( n ) (i-38b) 

The  blade-to-blade  periodicity  conditions  on  density  and  the  physical 
velocity  components  then  have  the  form  (cf.  Eq . (1-8)): 
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It  is  possible  to  construct  an  implicit  solution  scheme  for  the 
Navier-Stokes  equations  at  a cascade  periodic  boundary;  however,  such  a 
scheme  would  entail  rather  complicated  changes  to  the  Navier-Stokes  computer 
code.  Hence,  in  the  present  effort  an  explicit  scheme  has  been  adopted.  Thus 
the  form  of  the  governing  equations  applied  at  a periodic  boundary  follows 
from  Eq.  (1-16)  with  8=0;  i.e., 

Af"**  : At  ♦ S°  ] (I -AO) 


In  addition,  since  the  radial  coordinate  lines  of  the  cascade  coordinate  sys- 
tem do  not  possess  continuous  derivatives  at  periodic  boundaries  (see  Figs. 

10  and  11)  the  use  of  one  sided  difference  approximations  for  derivatives 
in  the  radial  coordinate  direction  is  indicated.  If  necessary,  an  implicit 
solution  scheme  using  central  differences  could  be  adopted  in  future  studies. 
Construction  of  an  implicit  scheme  would  involve  straightforward,  albeit 
very  tedious,  coding  changes.  The  introduction  of  central  difference 
approximations  in  the  radial  direction  is  contingent  on  future  refinements 
in  the  cascade  coordinate  system  described  in  Ref.  (1-16). 

With  the  foregoing  Ideas  in  mind  the  Navier-Stokes  computer  code  has 
been  modified  to  encompass  the  following  solution  scheme  for  the  blade 
passage  problem.  Intermediate  solutions,  ip  , of  the  direction  - 1 equations 
are  obtained  in  the  counterclockwise  direction  along  the  coordinate  loops 
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with  the  exception  of  the  blade  surface  and  outermost  loops.  These  solutions 
start  from  the  first  radial  coordinate  line  which  terminates  on  the  front 
end  cap  (i.e.,  at  point  A in  Fig.  11)  and  proceed  around  the  blade  to  the 
radial  coordinate  line  which  terminates  at  the  end  of  the  upper  periodic 
boundary  (i.e.,  point  H in  Fig.  11).  Solutions  at  the  n + 1 time  level, 

1 1>  , are  then  obtained  by  solving  the  direction  - 2 equations  along  the 

radial  coordinate  lines  subject  to  solid  surface  boundary  conditions  at 
the  blade,  upstream  and  downstream  far-field  conditions  on  the  front,  AB  in 
Fig.  11,  and  rear,  EF,  endcaps,  respectively,  of  the  outer  coordinate  loop 
ABCDEFGH  (Fig.  11),  explicit  Navier -Stokes  solutions  on  the  lower  periodic 
boundary,  CD,  and  blade-to-blade  periodicity  conditions  on  the  upper 
periodic  boundary  GH.  Solutions  along  radial  lines  which  terminate  at 
periodic  points  on  the  outer  boundary  proceed  simultaneously.  Although 
this  is  not  strictly  required  in  the  present  formulation,  which  uses 
explicit  Navier-Stokes  solutions  at  the  lower  periodic  boundary,  it  has 
been  incorporated  so  that  possible  future  improvements,  i.e.,  implicit 
procedures  at  the  lower  boundary  can  be  readily  incorporated. 

As  part  of  the  current  effort,  the  modifications  to  the  matrix  inversion, 
boundary  condition,  general  coefficient,  and  spatial  differentiation  sub- 
routines of  the  implicit  time-marching  Navier-Stokes  computer  code,  to 
incorporate  the  foregoing  formulation  have  been  completed.  The  changes  to 
the  existing  code  have  been  accomplished  in  a concise  manner  involving  the 
addition  of  only  one  new  subroutine  and  only  a limited  number  of  additional 
input  parameters.  The  code  modifications  have  been  checked  out  by  calculating 
the  symmetric  flow  past  an  isolated  circular  cylinder  using  symmetry  condi- 
tions, which  are  similar  in  form  to  cascade  periodicity  conditions  (i.e., 
symmetry  involves  only  a change  of  sign  in  Eq.  (I-39c))  at  the  endpoints  of 
radial  lines  extending  above  and  below  the  cylinder.  Computed  results  for 
this  case  after  five  time  steps  are  virtually  identical  to  those  obtained 
using  the  usual  isolated  body,  far-field  boundary  conditions  for  the  circular 
cylinder  calculations.  Such  agreement  indicates  that  the  modifications  to 
the  Navier-Stokes  code  required  to  treat  cascade  periodicity  conditions  have 
been  successfully  completed  and  that  this  code  is  ready  to  accept  the  cascade 
geometry  module  developed  by  Eiseman  (Ref.  1-16).  However,  before  the 
geometry  module  is  incorporated  into  the  Navier-Stokes  code  the  coordinate 
generator  must  be  further  developed  to  provide  the  higher  order  smoothness 
required  for  viscous  flow  studies. 


29 


REFERENCES  - PART  I 


1-1.  Briley,  W.  R.  and  H.  McDonald:  An  Implicit  Numerical  Method  for  The 
Multidimensional  Compressible  Navier-Stokes  Equations.  United 
Aircraft  Research  Laboratories  Report  M911363-6,  November  1973. 

1-2.  Douglas,  J.  and  J.  E.  Gunn:  A General  Formulation  of  Alternating 
Direction  Methods.  Numerische  Math.,  Vol.  6,  pp . 428-453,  1964. 

1-3.  Briley,  W.  R.,  H.  McDonald,  and  H.  J.  Gibeling:  Solution  of  the 

Multidimensional  Navier-Stokes  Equations  by  a Generalized  Implicit 
Method.  United  Technologies  Research  Center  Report  R75-911363-17 , 
January  1976. 

1-4.  Richtmeyer,  R.  D.  and  K.  W.  Morton:  Difference  Methods  for  Initial 
Value  Problems.  John  Wiley  & Sons,  1967. 

1-5.  Gibeling,  H.  J.,  S.  J.  Shamroth  and  P.  R.  Eiseman:  Analysis  of 
Strong-Interaction  Dynamic  Stall  for  Laminar  Flow  on  Airfoils. 

United  Technologies  Research  Center  Report  R77-912030-17 , July  1977. 

1-6.  Aris,  R.:  Vectors,  Tensors,  and  the  Basic  Equations  of  Fluid 
Mechanic s . Prentice-Hall,  Inc.,  1962. 

1-7.  Isaacson,  E.,  and  H.  B.  Keller:  Analysis  of  Numerical  Methods. 

John  Wiley  and  Sons,  Inc.,  New  York,  1966. 

1-8.  Roache,  P.  J.:  Computational  Fluid  Dynamics.  Hermosa  Publisher, 
Albuquerque,  New  Mexico,  1972. 

1-9.  Beam,  R.,  and  R.  F.  Warming:  An  Implicit  Finite  Difference 

Algorithm  for  Hyperbolic  Systems  in  Conservation-Law-Form.  Journal 
of  Computational  Physics,  Vol.  22,  pp.  87-110,  September  1976. 

1-10.  Coutanceau,  M. , and  R.  Bouard : Experimental  Determination  of  the 
Main  Features  of  the  Viscous  Flow  in  the  Wake  of  a Circular 
Cylinder  in  Uniform  Translation.  Part  1.  Steady  Flow.  Journal  of 
Fluid  Mechanics,  Vol.  79,  pp . 231-256,  1977. 

1-11.  Roberts,  G.  0.:  Computational  Meshes  for  Boundary  Layer  Problems. 
Proceedings  of  the  Second  International  Conference  on  Numerical 
Methods  in  Fluid  Dynamics,  Springer  Verlag,  New  York,  p.  171,  1971. 

1-12.  Kawaguti,  M. : Numerical  Solution  of  the  Navier-Stokes  Equations  for 
the  Flow  Around  a Circular  Cylinder  at  Reynolds  Number  40.  Journal 
of  the  Physical  Society  of  Japan,  Vol.  8,  pp.  747-757,  1953. 


30 


REFERENCES  - PART  I (Cont 'd) 


1-13.  Apelt,  C.  J.:  The  Steady  Flow  of  a Viscous  Fluid  Past  a Circular 
Cylinder  at  Reynolds  Numbers  40  and  44.  Aeronautical  Research 
Council  R.  & M.  No.  3175,  1958. 

1-14.  Kavaguti,  M. , and  P.  C.  Jain:  Numerical  Study  of  a Viscous  Fluid 
Flow  Past  a Circular  Cylinder.  Journal  of  the  Physical  Society 
of  Japan,  Vol . 21,  p.  2055. 

1-15.  Son,  J.  S.  and  T.  J.  Hanratty:  Numerical  Solution  for  the  Flow 

Around  a Cylinder  at  Reynolds  numbers  of  40,  200  and  500.  Journal 
of  Fluid  Mechanics,  Vol.  35,  pp . 369-386,  1969. 

1-16.  Eiseman,  P.  R. : A Coordinate  System  for  a Viscous  Transonic  Cascade 
Analysis.  United  Technologies  Research  Center  Report,  February  1977. 


31 


PART  II 


A VISCID/ INVISCID  INTERACTION  APPROACH  j 

FOR  HIGH  REYNOLDS  NUMBER  TRAILING 
EDGE  FLOW 

\ 

I 

. 1 

I 

1 


LIST  OF  SYMBOLS  - PART  II 


a 


C 


e 


h 

H 

H 

j 

3 

M 


n 

N 

£ 

P 

P 

q 


reference  length 

Chapman-Rubes  in  viscosity  constant 
constant  pressure  specific  heat 
unit  vector 

coordinate  scale  factor 
total  enthalpy 

surface  height  in  bottom  deck  variables 
geometry  index:  j = 0 2-D  flow;  j * 1 axisymmetric  flow 
Cartesian  coordinate  unit  vector 
Mach  number 

transposed  normal  coordinate 

transposed  stretched  normal  coordinate 

asymptotic  ordering  symbol 

static  pressure 

inner  deck  pressure 

heat  transfer 

transverse  coordinate  and  surface  radii  respectively 
Reynolds  number 

transposed  longitudinal  coordinate 

position  of  coordinate  system,  inner  variables 

position  of  coordinate  system 
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static  temperature 
longitudinal  velocity 

inner  deck  or  inviscid  longitudinal  velocity 

normal  velocity,  transposed  coordinates 

inner  deck  normal  velocity 

normal  velocity 

longitudinal  coordinate 

inner  deck  longitudinal  coordinate 

normal  coordinate 

inner  deck  normal  coordinate 

inner  region  ramp  angle  parameter 

ratio  of  specific  heats 

inner  deck  displacement  thickness 

perturbation  parameter 

longitudinal  surface  curvature 

Blasius  constant 

viscosity  coefficient 

inclination  angle  with  horizontal  axis 
density 

Prandtl  number 
shear  stress 
trailing  edge  angle 
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Subscripts 


e 

inviscid  edge  values 

F.P. 

flat  plate  solution 

m 

matching  region  values 

T.E. 

trailing  edge  location 

V 

viscous  or  main  deck  region 

W 

wall  values 

x.  y.  z 

coordinate  direction 

CD 

free  stream  value 
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GENERAL  APPROACH 


The  problem  of  base  flow  separation  off  the  trailing  edge  region  of 
airfoils  (i.e.,  trailing  edge  stall)  causes  difficulty  in  analytically 
predicting  the  losses  and  loading  in  turbomachinery  appl icat ions . Upon  close 
inspection,  it  is  obvious  that  this  flow  and  its  attendent  difficulties  are 
very  similar  to  those  encountered  in  predicting  surface  separation  bubbles, 
the  main  difference  here  being  that  reattachment  takes  place  on  a slip  line 
instead  of  on  a solid  surface.  While  it  has  long  been  clear  that  classical 
boundary  layer  theory  is  inadequate  for  predicting  separated  flows,  recent 
analytical  developments  indicate  that  a consistent  and  reliable  analytical 
model  can  be  constructed  using  somewhat  less  than  the  full  Navier-Stokes 
equations.  Additionally,  careful  analysis  of  experimental  data  suggests  that 
the  effect  of  the  stalled  flow  trailing  edge  separation  bubble  on  the  flow 
field  is  confined  to  the  immediate  region  in  and  about  the  trailing  edge 
region.  This  observation  indicates  that  to  capture  the  trailing  edge  flow 
phenomenon  it  is  only  necessary  to  modify  the  classical  approach  (i.e., 
inviscid  flow  plus  boundary  layer  theory)  in  the  immediate  vicinity  of  the 
trailing  edge.  Thus,  while  methods  based  on  the  application  of  the  full 
Navier-Stokes  equations  to  the  entire  flow  field  will  surely  be  formally 
applicable,  alternate,  simpler  approaches  are  now  available. 

For  separation  bubble  flows,  all  of  the  simpler  analytical  models  to 
date  contain  two  common  features:  first,  in  the  viscous  region  diffusion 
effects  normal  to  the  streamlines  dominate  all  others  and;  second  the  inviscid 
flow  field  represents  flow  over  a surface  formed  by  thickening  the  original 
shape  with  the  local  viscous  region's  displacement  thickness.  Variations  on 
this  approach  involve  either  a composite  of  these  two  layers  (the  thin  layer 
approximation,  the  parabolized  approach,  etc.)  or  a substructure  delineation 
(the  triple  deck,  or  the  multi-deck  approaches).  In  either  case,  successful 
modeling  of  separation  bubble  flows  has  been  achieved  and  it  remains  now  to 
apply  and  extend  these  concepts  to  stalled  subsonic  and  transonic  trailing 
edge  regions  in  t urbomach inery  appl icat ions . This  is  the  overall  goal  of  the 
present  study  with  specific  interest  in  the  formulation  of  the  two  layer, 
interacting  laminar  boundary  layer  model  of  this  flow  field.  A companion 
study  has  also  been  conducted  on  the  development  and  demonstration  of  numer- 
ical techniques  for  solving  such  flows  for  the  supersonic  trailing  edge 
problem  - the  results  of  which  will  be  summarized  here. 

The  general  approach  to  be  employed  here  follows  closely  the  work 
of  Melnik,  Chow,  and  Mead  (Ref.  I I— 1 ) who  developed  an  interacting  boundary 
layer  model  for  high  Reynolds  number  flow  past  isolated  airfoils.  As  depic- 
ted in  Fig.  12,  the  overall  flow  field  model  is  rather  conventional  except 
that  an  interaction  model  is  employed  to  bridge  the  trailing  edge  region 
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where  the  blade  boundary  layer  transits  to  a wake  flow.  The  major  element 
of  this  approach  is  that  for  high  Reynolds  number  the  laminar  viscous  flow  is 
assumed  to  influence  the  inviacid  flow  field  principally  through  displacement 
effects  and  is  represented  by  the  displacement  thickness  as  added  to  the  air- 
foil surface  and  wake  centerlines.  Thus,  the  major  elements  of  the  flow  are: 
(1)  an  inviscid  flow  outside  the  thin  viscous  regions  here  taken  to  be  gov- 
erned by  the  potential  flow  equations  as  applied  to  cascades;  (2)  a weakly 
interacting  boundary  layer  region  over  the  forward  portion  of  the  blade;  (3) 
a weakly  interacting  wake  region  aft  of  the  trailing  edge  region;  and  (4)  a 
strongly  interacting  trailing  edge  region  where  the  local  pressure  levels 
must  be  established  along  with  the  viscous  solution.  Coupling  of  these  four 
flows  produces  a boundary  value  problem  in  the  flow  direction  that  Melnik 
et  al . (Ref.  11—1)  have  shown  can  be  solved  with  iterative  techniques. 
Attention  here  will  be  directed  at  the  solution  of  the  strongly  interacting 
trailing  edge  region  where  Melnik  et  al.  used  approximate  techniques  to 
solve  for  attached  flows  only.  The  separated  (or  stalled)  trailing  edge 
problem  will  be  addressed  here  using  recently  developed  interacting  boundary 
layer  concepts  (see  Refs.  11-2,  6).  Attention  will  be  focused  on  the  laminar 
flow  pa6t  symmetric  6harp  trailing  edges,  as  depicted  in  Fig.  13,  with  trail- 
ing edge  angles  large  enough  to  cause  the  boundary  layer  to  separate  off  the 
surface  and  reattach  on  a wake  slip  line. 
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THE  TRAILING  EDGE  FLOW  MODEL 


The  details  of  the  flow  field  for  the  trailing  edge  region  are  given  in 
Fig.  13.  Here  the  basic  inviscid  flow  is  that  of  flow  into  and  out  of  a 
concave  corner.  This  strongly  deaccelerat ing  flow  field  forces  the  incoming 
boundary  layer  to  grow,  deflect  the  inviscid  stream,  and  finally  separate 
from  the  surface  and  reattach  along  the  wake  - producing  a displacement 
surface  that  smooths  the  sharp  corner's  influence  on  the  local  inviscid  flow 
field  ( see  Fig . 14) . 

The  analytical  model  for  this  flow  field  is  taken  directly  from  that 
already  proved  to  be  applicable  for  flow  into  a compression  corner  made  up  of 
two  6olid  walls.  Vatsa  and  Werle  (Ref.  11-4),  Jenson,  Burggraf,  and  Rizzetta 
(Ref.  1 1— 7 ) and  Burggraf  et  al  (Ref.  1 1— 11)  have  clearly  shown  that  the 
interacting  boundary  layer  concept  provides  a rational  approximation  for  such 
flows  and  this  concept  will  be  carried  over  directly  here.  The  principle 
governing  equations  for  the  viscous  region  are  formally  recovered  from  the 
Navier-Stokes  equations  using  the  formalism  of  higher  order  boundary  layer 
theory  (Ref.  I I— 8 ) to  identify  the  displacement  thickness  corrections  to  the 
flow  equat ions . 

The  basic  approach  taken  here  is  similar  to  that  first  set  down  by  Van 
Dyke  (Ref.  1 1— 9 ) except  that  here  no  attempt  is  made  to  filter  out  the 
higher-order  effects  into  separate  linear  problems.  Van  Dyke  has  given  a 
statement  of  the  full  Navier-Stokes  equations  for  two-dimensional  or  axisym- 
metric  flow  in  terms  of  the  coordinates  x and  y with  all  distances  referenced 
to  a length,  a,  velocities  to  , pressure  to  p.U^2,  density  to  p^, , temperature 
to  U^/Cp,  enthalpy  to  U.2,  and  viscosities  co  the  value  of  u at  T ■ U.VCp. 
With  these  definitions  the  characteristic  Reynolds  number  becomes 


Re  - P.Ua,a/p(U-2/Cp)  - 1/e8 


(II-l) 


where  e becomes  the  principle  perturbation  parameter  for  an  asymptotic 
analysis  of  the  separated  flow  region  using  the  concepts  developed  by 
Stewartson  (Ref.  11-10). 

For  present  purposes  we  approach  the  full  equations  with  the  assumption 
that  in  a vanishingly  small  region  near  the  surface,  a region  of  order 
thick,  all  flow  properties  except  the  normal  velocity  component  v are  of 
order  one  - v itself  being  of  order  e^. 

Without  loss  of  generality  the  coordinate  system  used  here  will  be  placed 
placed  along  a smooth  surface  always  located  near  the  line  of  zero  longitudinal 
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velocity  as  shown  in  Fig.  15.  It  will  be  shown  later  that  within  the 
framework  of  an  analytical  model  applicable  for  asymptotically  large 
Reynolds  number,  this  arbitrariness  in  the  coordinate  definition  will 
i be  totally  acceptable. 

The  governing  equations  are  recovered  by  introducing  a stretching  of  the 
viscous  region  normal  coordinate,  y v , to  give  a boundary  layer  coordinate  y 
as 

y = yj//«  ( i i-2a  > 


and  retaining  the  longitudinal  coordinate  scale  as 


x = X„ 


( 1 1—  2b ) 


and  using  unsubsc r ipted  variables  to  designate  the  flow  properties  in  the 
viscous  region  as  follows: 


longitudinal  velocity: 

u ytXj/,  y € ) = u(x,y*€) 

( I I-3a) 

normal  velocity: 

(II-3b) 

dens ity : 

^xi/.yj/i«)  =A>Cx,yi«) 

(II-3c) 

pressure : 

Pi/(xi/.yi/i<)=p  (*.y.«) 

( I I-3d ) 

static  enthalpy: 

sTU,y,«) 

( I I-3e ) 

total  enthalpy: 

Hi/lxi/>y„i«)  = H(x,yi«) 

( 1 1—3  f ) 

The  coordinate  scale  factors 

are  designated  as 

4 

hx=  | + <4*y  sh 

( I I-4a ) 

hy  = l 

(II-Ab) 

i 

h 

z = (r0  + c4y  cos#)*  * r1 

( I I-4c ) 

Introducing  these  definitions  into  the  full  equations,  keeping  all  terms  up 
to  second  order,  and  invoking  the  asymptotic  matching  principle  to  determine 
the  interaction  with  the  inviacid  flow  field,  gives  the  following  for  plane 
(j  ■ 0)  or  axisymmetric  flows  (j  » 1). 
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continuity : 


longitudinal  momentum: 


4r 

ax  [Jy 


where 


total  enthalpy: 


P 


_d_  (pr^u)  + (phr J w)  = 0 (Il-5a) 

ax  aN 


au 

ax 


ahu 

ay 


-?m(u«/h)a7(Ue/h)+ 


ph2u2)dy 


_L  X (rih2r) 

hr)  a y 


(II-5b) 


”'*(77 -l,h»/h) 


(II-5c) 


JL 

dy 


hr  • (q  + ur) 


(II-5d) 
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where 


q s/x/cr 


dT 

77  = /x/a 


(3H  t3  u 

d y d y 


( I I-5e ) 


Also  the  viscosity  law  is  given  as: 


M = f(T) 


where  f(T)  represents  any  appropriate  temperature  viscosity  law  with 


T = H - UZ/2 


( 1 1-6 ) 


The  final  two  relations  needed  are  derived  from  the  matching  conditions  as 


p-  +/  (pmu«J -/>"'u*]dy  y/>(h-t) 


(Il-7a) 


and 


Pm 


,JL 

/pe  = (/5m/pe)y  = [l+(Ue2/2)(He-Ue2/2)(l-l/h2)Jy-1  (II-7b) 


The  edge  conditions,  Ue  and  Hg  identified  above,  are  to  be  obtained 
from  the  inviscid  flow  field  corrected  for  displacement  thickness  effects. 
This  issue  will  be  addressed  later  in  this  section  with  attention  here 
directed  toward  application  of  the  above  general  equations  to  the  trailing 
edge  problem  and  identification  of  appropriate  boundary  conditions. 
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The  first  simplification  invoked  here  for  the  trailing  edge  problem 
involves  restriction  to  the  two  dimensional  case,  i.e.,  taking  j = 0 in  the 
above  equations.  The  second  simplification  comes  about  through  elimination 
of  the  curvature  induced  effects  identified  through  the  scale  factor  h. 

This  is  accomplished  on  two  bases.  First,  in  regions  of  weak  interaction, 
such  as  far  ahead  and  aft  of  the  trailing  edge,  the  geometry  of  the  coordi- 
nate line  is  flat  and  thus  its  curvature  is  zero.  Second,  in  regions  of 
strong  interactions  as  occur  near  the  trailing  edge  point,  the  formal  ordering 
analysis  for  asymptotically  large  Reynolds  number  indicates  that  these  terms 
will  be  of  secondary  importance.  To  demonstrate  this  point  it  is  first  neces- 
sary to  show  that  the  strong  interaction  equations  derived  from  the  asymptotic 
triple  deck  approach  are  subsets  of  the  interacting  boundary  layer  equations 
presented  above.  Burggraf,  Rizzeta,  Werle  and  Vatsa  (Ref.  11-11)  have  proven 
this  point  numerically,  while  Vatsa  (Ref.  11-12)  and  Vatsa  and  Werle  (Ref. 

11-4)  have  shown  that  the  introduction  of  the  triple  deck  scaling  laws  into 
the  interacting  boundary  layer  equations  leads  directly  to  the  correct  subscale 
equations  for  strong  interaction  regions. 

Thus,  for  example,  for  trailing  edge  angles  a “(9(c^)  it  is  found  that 
within  a region  x - xTE  “ , y =$(e)'.  the  bottcxn  deck  region 

within  which  the  separation  bubble  is  contained,  the  dependent  variables 
scale  as 


u = 0 (e) 

( I 1-8 a) 

w =0(l/«) 

( I I— 8b ) 

Ap  =<%  (€2) 

( I I-8c ) 

p = <%(  1) 

►— « 

1 

00 

CL 

H=0(!) 

( I I-8e ) 

fJL.  =(*0) 

( 1 1— 8 f ) 

In  addition,  since  in  the  present  approach  the  coordinate  axis  is  placed  on  a 
surface  necessarily  contained  within  this  region,  it  follows  that,  the 
coordinate  location,  tQ(x)  ■#(e'*)  and  its  curvature  r »$(l/e).  With 
these  scalings,  the  dominate  terms  of  the  governing  equations  for  the  bottom 
deck  can  be  identified  for  the  continuity  equations  as 


cont inui t y : 


( I l-9a) 


See  the  following  section  of  this  report  for  a detail  description  of  the 
scaling  laws  discussed  here. 
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energy : 


P. 

a 


( I I-9b ) 


u±i 

d x 


+ w 


dH 

dy 


The  limit  form  of  the  longitudinal  momentum  equation  is  dependent  on  the  scaling 
law  for  the  edge  velocity  Ue , which  is  established  in  the  outer  most 
portion  of  the  strong  interaction  region  (y  “£(e^)).  in  this  region, 
one  finds  that  the  viscous  perturbation  to  With  this  the 

longitudinal  momentum  Eq.  II-5b  reduces  to 


2 

du  du  dUe  d u 

— — + w — Ue  — — - 

d x d y d x dy  z 


( II-9c ) 


A similar  study  of  the  appropriate  scalings  in  the  middle  deck  region  of 
the  strong  interaction  process  would  again  reveal  that  the  curvature  terms 
were  of  secondary  importance,  but  that  here  density  and  temperature  variations 
were  of  lead  order.  Thus  a composite  set  of  equations  can  be  established  for 
the  two  regions  as 


+ =0 
d x d y 


( II- 10a ) 


’ujLy_ 

du 

+ w 

- p u dU*  = 

d 

L 

(II-10b) 

d x 

dy 

e dx 

dy 

v dy  I 

, dH  __  dH 

. _d_ 

P dT  ^ „„  du 

P 

u—  + w— 
dx  dy 

dy 

dy  hdy 

(ll-10c) 


which  are  merely  the  classical  boundary  layer  equations  with  edge  properties 
adjusted  for  interaction  effects. 

The  boundary  conditions  for  the  present  equations  can  be  set  down 
directly  as 
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y — ► ao  ; u— ^ Ue  and  H — ► He 

(11-11) 

y = - t0(x)/«4  * - t(x) 

( 1 1- 1 2a ) 

On  the 

blade  surface:  x £ XT  E . ’ u “ « * 0.  H “ Hw 

( I I- 1 2b ) 

On  the 

• v i • V = 0 

wake  centerline:  x > xT  F ; 1 

( 1 1-  12c  ) 

' - j • 7H  = 0 

( 1 1 - 1 2d  ) 

J-?(7v)  = 0 

( 1 1- 1 2e ) 

where  j is  the  vertical  unit  vector. 


Following  Jenson  et  al  (Ref.  11—7)  these  latter  conditions  can  be 
considerably  simplified  through  use  of  Prandtl's  Transposition  Theorem.  Thus 
de  f ine 


N = y ♦ t(x)  = n/c4 


(1 1-1 3a) 


so  that  the  innermost  boundary  conditions  are  placed  at  N * 0 * n. 
In  addition  define 


S = X 


( 1 1 - 1 3b) 


and 


v = w+u^-  = w-f  ut' 
ax 


so  that 


j_ _ _a  . t/jL 
« * » . * * 


<3  x <3s 


<JN 


( 1 1- 1 3c  ) 


( 1 1- 14a  ) 


and 


(II-Kb) 

dy  dN 


44 


With  these  definitions  the  governing  equations  remain  unchanged  in  form  and 
thu6  become 


dpu  dpv 

as  * ~dH 


( 1 1-1 5a  ) 


and 


P 


+ V 


Al 

du 


( 1 1-1 5b ) 


X 

X 

. _a_ 

p di  ..  au 

p 

u as  v aN 

‘ aN 

<y  aN  " ^ aN 

( 1 1 — 1 5c  ) 


The  boundary  conditions  for  these  equations  are  straightforward  on  the 
surface,  but  can  be  considerably  simplified  over  those  given  in  Eqs.  (11  — 1 2b ) 
through  (11—1 2d ) for  the  wake  centerline. 

Referring  to  Fig.  15,  consider  the  condition  that 


j-^  = 0 


(II— 16a) 


or 


w*sin<£-  u„cos<£  = 0 


(1 1-1 6b  ) 


Thus 


,=  «4  w 


,ton<£=  - ~ = -«Vu 


( 1 1- 1 6c  ) 


which  from  Eq.  ( 1 1 - 1 3c  ) give  that  along  the  alip  line, 


fl  = 0 V =0 


(11  — 1 6d ) 
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Evaluation  of  the  normal  derivative  conditions  is  somewhat  more  tedious. 
The  gradient  operator  is  given  as 


y 


h 


(II— 17a) 


so  that 


-y  : rang. 
1 h 


af ' '°"e  £ 


(11-1 7b) 


Thus  the  wake  centerline  condition  on  total  enthoipy  becomes 


dH  _ _ t0  dH 

dn  h + 10'?  ds 


( 1 1- 1 8a ) 


and  since 


7-  V = u/COS  8 


( 1 1-1 8b ) 


then 


dll  . 
dn 


tp 

h+t0'2 


<3u 

ds 


-IttftU 


( 1 1-1 8c ) 


These  relations  can  be  simplified  by  noting  that  in  the  weak  interaction 
region  tQ  -►  0 whereas  in  the  strong  interaction  region,  the  bottom  deck 
scaling  laws  apply.  Applying  Eqs.  (ll-8a)  through  (Il-8f)  to  Eqs.  (Il-17a) 
and  ( 1 1—  17b)  leads  to  the  result  that  for  o * 0,  s > s-pj-.  the  appropriate 
boundary  conditions  are  given  as 
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(11-19) 


dH 

an 


= 0 


du 

an 


= 0 


(II-20) 


It  remains  now  only  to  establish  the  analytical  basis  for  establishing  the 
edge  properties,  Ug,  pg , etc.  Both  the  weak  and  strong  interaction 
theories  indicate  that  this  can  be  accomplished  employing  the  effective 
displacement  surface  (see  Refs.  11-4,12,13,  and  15  for  example)  obtained 
by  the  simple  addition  of  the  boundary  layer  and/or  wake  displacement  thickness 
to  the  original  inviscid  flow  shape.  This  is  the  approach  that  will  be 
employed  here . 
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THE  HIGH  REYNOLDS  NUMBER  PROBLEM  - THE  TRIPLE  DECK  EQUATIONS* 


While  the  interacting  boundary  layer  problem  formulated  above  is  expected 
to  have  application  over  a very  wide  range  of  Reynolds  numbers,  it  has  been 
shown  by  Vatsa  and  Werle**  (Ref.  II-4)  that  these  equations  reduce  to  a rela 
tively  simple  and  instructive  form  as  the  Reynolds  number  increases.  In 
particular,  Vatsa  and  Werle  were  able  to  show  that  the  interacting  boundary 
layer  equations  contain  the  subscale  structure  (the  triple  deck  structure) 
identified  by  Stewartson  (Ref.  11-10),  and  Jenson,  Burggraf  and  Rizzeta  (Ref. 
II-7)  for  interacting  flows.  Of  particular  interest  here  is  the  application 
of  this  concept  to  the  supersonic  compression  corner  problem  as  depicted  in 
Fig.  16.  Through  identification  of  the  subscalings  shown  near  the  hinge  line 
and  of  the  ramp  angle,  a distinguished  limit  solution  to  the  Navier-Stokes 
equations  for  Re  >>  1 was  identified  that  captured  the  interaction  and  separa- 
tion phenomena.  As  such  the  solutions  are  as  fundamentally  important  to  the 
field  of  separation  flow  theory  as  the  Blasius  solution  is  to  attached  flow 
theory.  An  additional  benefit  of  this  approach  is  that  it  is  found  to  reduce 
the  parametric  dependence  of  the  problem  from  four  (M  , Re^,  T /T^,  6) 
to  one,  a,  where 


-1/2 


C(M„Z-|) 


1/4 


(11-21) 


Thus  these  equations  provide  a convenient  and  useful  test  vehicle  for  develop- 
ing solution  methods  and  demonstrating  concepts  while  providing  meaningful 
quantitative  separated  or  stalled  flow  solutions. 

The  approach  used  by  Vatsa  and  Werle  (Ref.  II-4)  follows  directly  frcro 
that  of  Jenson  et  al  (Ref.  II-7)  in  that  the  scaling  laws  identified  for 
strong  interactions  by  Stewartson  (Ref.  11-10)  are  applied  in  the  near 
region  of  a geometric  obstruction  - in  this  case,  a compression  corner.  The 
principle  region  of  interest  here  is  the  bottom  deck  of  Fig.  16,  which  can  be 
identified  in  terms  of  scaled  dependent  variables,  X and  Y,  given  as: 


x-xTE=  s-sTE 


X 


( I I-22a ) 


The  work  presented  here  and  in  the  following  two  sections  concerning  the 
asymptotic  Triple  Deck  equations,  their  numerical  solution  algorithm,  and 
the  demonstration  of  results  were  performed  under  a United  Technologies 
sponsored  project  outside  the  scope  of  the  contract  requirements.  They  are 
included  here  to  demonstrate  the  utility  of  this  approach. 

These  results  have  also  been  discussed  by  Burggraf  (Ref.  II-6)  and  will  be 
appearing  shortly  in  Ref.  II— 1 1 . 
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( 1 1 — 2 2b  ) 


n = y*+  <0  = €<N 


X-3/4 


[l 9 l/e] — (^/t,,)372  y 

(M*2-!)'/8 


The  dependent  variable  scalings  for  this  region  have  been  discussed  briefly 
in  an  earlier  section  of  this  report  and  are  given  here  in  detail  as: 

longitudinal  velocity: 


= \ 1/4 


U*  = U * X 


€ C 


1/8 


(M„  - I) 


1/8 


(t„/tw)i/2u 


normal  velocity: 


W^  + U^to'  = €4v 


(Tw/T«)l/2  V 


( 1 1-2  3a  ) 


( I l-23b  ) 


static  pressure: 


P 


P®  * 


€ C 


1/8 


( I I-23c ) 


With  the  introduction  of  these  asymptotic  scale  laws  and  those  for  the  local 
inviscid  flow  into  the  interacting  boundary  layer  equations  given  by  Eqs. 

( 1 1-1 5a ) through  ( 1 1- 1 5c  ) above,  Vatsa  and  Werle  (Ref.  11-4)  were  able  to 
identify  the  dominant  terms  of  the  interacting  boundary  layer  equations  for 
high  Reynolds  number  strong  interactions  as 

cont inuity : 


du 

d* 


(II-24a) 
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longitudinal  momentum: 


+ V 


au 

dx 


+ 


dL  = a!u 

dy  d'f 2 


( 1 3 -24b ) 


and  the  energy  equation  is  found  to  be  uncoupled  from  these  fundamental 
equations.  These  equations  are  identical  to  those  derived  from  the  Navier- 
Stokes  equations  by  Jenson  et  al  (Ref.  II-7).  The  viscous/invi sc  id  interac- 
tion law  for  the  asymptotic  problem  is  a simplified  version  of  that  for  the 
interacting  boundary  layer  problem  and  is  given  as 


P 


(II-25a) 


where  Hg  is  the  surface  height  given  here  for  the  ramp  in  inner  deck 
variables  as 


Hs=0  for  X $ 0 


( II  — 25b) 


H = a x for  X > 0 


( I l-25c ) 


and  s i6  a displacement  thickness  like  function  given  as 


8 = lim  (y-u) 

Y—  CD 


(II-25d) 


It  has  been  demonstrated  that  these  equations  provide  the  same  solution 
at  high  Reynolds  number  as  that  of  the  interacting  boundary  layer  equations 
by  Burggraf,  Rizzeta,  Werle  and  Vatsa  (Ref.  1 1 — 1 1 ) with  the  significant 
results  repeated  here  for  completeness.  For  a fixed  value  of  a * 2.5,  Figs. 
17  and  18  compare  separated  flow  solutions  of  the  interacting  boundary  layer 
equations  with  solutions  of  the  fundamental  equations  given  above. 


Also  see  Refs. 


11-4,  6,  and  11 


for  further  discussion  of  this  point. 
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It  is  clearly  seen  that  both  the  surface  pressure  and  skin  friction 
distributions  obtained  from  the  interacting  boundary  layer  approach  system- 
atically approach  the  triple  deck  solutions  as  the  Reynolds  number  increases  - 
thus  verifying  the  utility  of  the  fundamental  triple  deck  version  of  the 
governing  equations  for  study  of  the  interacting  boundary  layer  model  of 
the  trailing  edge  stall  problem.  For  this  case  the  boundary  condition  aft  of 
the  trailing  edge  is  obtained  by  applying  the  scaling  laws  above  to  Eqs. 

( 1 I— 1 6d ) and  ( II— 18d)  to  obtain: 


X > 0,  Y = 0 


V = 0 


( II-26a) 


M --o 

av 


( I l-26b ) 


To  provide  insight  into  the  nature  of  the  trailing  edge  flow  field,  solutions 
of  Eqs.  (11-24)  through  (11-26)  have  been  obtained  using  numerical  techniques 
outlined  in  the  following  section  with  a discussion  of  the  results  presented  in 
the  subsequent  section. 
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NUMERICAL  METHOD  - THE  TRIPLE  DECK  EQUATIONS 


The  numerical  method  applied  to  solve  the  trailing  edge  problem  described 
above  will  be  outlined  briefly  below.  The  approach  represents  an  extension 
of  the  concept  originally  developed  by  Werle  and  Vatsa  (Ref.  1 1—  3 ) for  the 
full  interacting  boundary  layer  equations.  The  method  used  here  directly 
follows  the  application  and  generalization  of  this  technique  to  the  triple 
deck  equations  by  Napolitano,  Werle,  and  Davis  (Ref.  1 1— 2 ) for  subsonic  and 
supersonic  separated  flow  past  protuberances  and  compression  corners.  The 
basic  approach  uses  a marching  type  implicit  finite  difference  scheme  to 
solve  the  continuity  and  momentum  Eqs.  (II-4a)  and  (II-4b)  with  nonlinear- 
ities approximated  using  previous  station  data.  The  pressure  interaction 
effect  enters  through  the  pressure  gradient  term  of  Eq.  (II-4b)  and  is 
coupled  to  the  local  velocity  through  Eqs.  (II-5a)  through  (II-5b).  This 
effect  is  accommodated  through  introduction  of  a time-like  relaxation  scheme 
wherein  dP/dX  of  Eq.  (II-4b)  is  replaced  with  the  term  dP/dX  - d6/dt  and  the 
solution  marched  in  time  to  achieve  d6/dt  =0.  A simple  superposition 
technique  is  employed  to  determine  the  simultaneous  solution  of  the  conserva- 
tion and  interaction  equations  at  a given  station  so  that  computational 
efficiency  is  maintained. 

For  completeness  two  versions  of  the  numerical  algorithm  were  employed 
here.  The  first  of  these  was  a slight  modification  of  the  numerical  algorithm 
presented  in  Ref.  1 1—  2 and  as  such  was  first  order  accurate  in  the  longitu- 
dinal direction,  AX.  A second  order  accurate  version  of  the  numerical  method 
was  also  available  and  has  been  applied  here.  A brief  outline  of  the  prin- 
ciple features  of  each  algorithm  is  given  below. 


The  First  Order  Method 


The  algorithm  employed  here  is  virtually  identical  to  that  developed 
by  Napolitano,  Werle,  and  Davis  (Refs.  1 1— 2 and  11-15)  for  subsonic  and 
supersonic  flows  with  several  changes  as  required  to  accommodate  the  de- 
tails of  the  flow  structure  in  the  trailing  edge  region.  These  changes  for 
the  present  supersonic  case  represent  generalizations  of  those  concepts 
employed  by  Napolitano  et  al , (Refs.  1 1—2  and  11-15)  to  accommodate  sharp 
corner  effects  in  subsonic  flows.  The  changes  are  identified  below  with 
reference  to  the  appropriate  sections  and  equations  of  Reference  11-15. 

(i)  the  continuity  equation  (II-24a)  solution  technique  was  modified 
one  station  aft  of  the  trailing  edge  in  accordance  with  Section 
IV- 2 and  Eq . (4.13)  of  Ref.  11-15. 
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(ii) 


the  numerical  representation  of 
(ll-24b)  was  modified  according 
Ref.  11-15. 


the  normal  convection  term  of  Eq . 
to  Section  1V-2  and  Eq . (4.14)  of 


(iii) 


at 

a 

given  time 

in 

Eq 

. ( I I-24b) 

t i 

ons 

us i ng  Eq  . 

level,  the  pressure  gradient  term  dP/dx  appearing 
above  was  obtained  from  previous  time  level  solu- 
(4.1)  of  Ref.  11-15. 


(iv)  the  numerical  approach  employed  in  Section  IV-4  of  Ref.  11-15  to 
solve  the  interaction  law  and  thus  update  the  pressure  gradient 
estimate  was  employed  here.  However,  the  Cauchy  integral  terms 
for  subsonic  flow  were  replaced  by  the  appropriate  algebraic  rela- 
tion for  supersonic  flows. 


A final  modification  was  made  to  provide  more  accurate  representation  of 
the  local  pressure  levels  after  the  solutions  were  completed.  This  was  found 
necessary  principally  for  the  flat  plate  trailing  edge  problem  where  singular 
gradients  and  thus  large  truncation  errors  were  encountered.  Thus,  the  pres- 
sure at  a typical  point,  point  2,  was  given  in  terms  of  previous  station 
values,  point  1,  from  Eq . ( II— 25a)  as 


(S2+h52)-(S,-hS|) 

AX 


(11-27) 


The  Second  Order  Method 


In  an  effort  to  provide  very  accurate  numerical  baseline  solutions  so 
that  a valid  independent  assessment  of  the  analytical  model  could  be  achieved 
untainted  by  numerical  complications,  an  available  second  order  accurate  (in 
AX)  numerical  algorithm  was  also  applied  to  the  present  problem.  The  algo- 
rithm employs  an  approach  very  similar  to  that  presented  in  Ref.  11-15  except 
that  here  the  first  order  accurate  implicit  windward  differencing  used  in  the 
longitudinal  direction  for  the  momentum  equation  was  replaced  by  a second 
order  accurate  Crank-Nicolson  scheme.  Similarly  the  accuracy  of  the  conti- 
nuity equation  solution  was  upgraded  to  second  order  in  AX  by  centering  the 
difference  equation  between  two  longitudinal  stations.  To  avoid  large  numeri- 
cal errors  due  to  anticipated  iump  discontinuities  in  the  flow  variables 

central  difference  approxima- 
a single  station  immediately  aft 

ot  the  trailing  edge. 


several  of  the  formally  second  order  accurate 
tions  were  relieved  to  a first  order  level  at 
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RESULTS  AND  DISCUSSION  - THE  TRIPLE  DECK  PROBLEM 


Solution  to  the  high  Reynolds  number  form  of  the  interacting  boundary  r 

layer  equations  (the  triple  deck  equations)  were  first  obtained  for  super- 
sonic flow  into  a compression  corner  as  depicted  in  Fig.  19.  As  such,  this 
geometry  only  differs  from  that  of  the  trailing  edge  in  that  the  center  line 
symmetry  condition  is  replaced  by  a zero  slip  condition  aft  of  the  corner 
point.  An  essentially  exact  solution  is  available  for  this  problem  from  the 
work  of  Rizzetta  (Ref.  11-13)*,  thus  allowing  a basis  for  assessing  the  cur- 
rent algorithms.  Figure  20  gives  a comparison  of  the  normalized  wall  shear 
stress  and  surface  pressures  for  a reduced  ramp  angle  of  a * 2.5  — a case 
for  which  separation  and  flow  recirculation  was  encountered.  As  shown  in  the 
figure  both  of  the  present  algorithms  provide  reasonable  approximations  to  the 
exact  results  showing  a pressure  rise  ahead  of  the  corner  up  to  a "plateau" 

level  as  the  flow  separates  over  the  corner  followed  by  a second  pressure  rise 

to  the  ramp  pressure  level  during  flow  reat tachment . Quite  naturally  it  is 
found  that  the  second  order  accurate  (in  AX)  algorithm  is  more  accurate  than 
its  first  order  accurate  counterpart.  These  solutions  were  obtained  on  a 
UNIVAC  1110  computer  in  approximately  30  seconds  of  CPU  time  for  grid  sizes  of 
70  (in  X)  by  26  (in  Y)  and  a convergence  criterion  that  the  average  variation 
per  iteration  be  less  than  10-l!*.  As  shown  in  the  figure  the  accuracy  was 
found  to  be  somewhat  insensitive  to  the  normal  grid  size  (AY)  apparently  due 

to  the  second  order  level  of  the  algorithm  for  the  normal  direction.  In  order 

to  verify  the  utility  of  the  present  algorithms,  solutions  were  also  obtained 
across  the  full  range  of  a with  the  resulting  surface  pressure  and  shear  dis- 
tributions given  in  Figs.  21  and  22  respectively.  These  results  clearly 
show  the  significant  and  systematic  growth  of  the  interaction  region  ahead  of 
the  corner  as  the  ramp  angle  is  increased. 

The  algorithm  as  employed  to  perform  the  ramp  calculations  above,  was 
modified  for  the  trailing  edge  problem  through  a single  adjustment  aft  of 
the  corner  point  to  replace  the  no  surface  slip  condition  with  the  wake  center- 
line  symmetry  condition  of  Eqs.  (11-26).  Solutions  were  then  obtained  for 
a range  of  a Up  to  and  including  the  stall  (separation)  condition.  Compari- 
sons are  first  given  here  for  the  flat  plate  (<*  “ 0)  case  since  Daniels  (Ref. 

11-14)  has  already  provided  very  accurate  solutions  for  this  case  which  can 
be  used  to  assess  the  accuracy  of  the  present  algorithms.  Figure  23  gives 
a detailed  assessment  of  the  present  results  obtained  from  both  the  first  and 
second  order  (in  AX)  algorithms.  Here  the  pressure  distribution  on  the  plate 
surface  up  to  the  trailing  edge  is  given  along  with  the  wake  centerline  pres- 
sure distribution.  For  this  geometry,  interaction  effects  cause  the  pressure 

See  also  Refs.  II-4,  6 and  11 
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to  begin  to  drop  far  ahead  of  the  trailing  edge  as  it  anticipates  the  abrupt 
, boundary  layer  thinning  aft  of  the  trailing  edge  when  the  retarding  influence 

of  the  surface  no  slip  condition  is  relieved.  Both  the  first  and  second  or- 
i der  results  essentially  reproduce  Daniels  (Ref.  11-14)  results  with  the  only 

significant  difference  occurring  in  the  immediate  vicinity  of  the  trailing 
edge  itself.  The  small  oscillations  observed  in  the  centerline  pressures  for 
the  second  order  algorithm  and  "kink"  observed  for  the  first  order  algorithm 
are  caused  by  the  extremely  high  acceleration  of  the  flow  immediately  aft  of 
the  trailing  edges.  As  shown  in  Fig.  24  these  minor  difficulties  are  elimi- 
nated as  the  grid  spacing  is  refined  with  the  resulting  1st  order  solution 
providing  a "smoothed"  version  of  the  second  order  results  near  the  trailing 
edge  . 

Solutions  for  the  general  trailing  edge  problem  were  then  generated 
l using  both  of  the  present  algorithms.  Figures  25  through  27  give  the  second 

order  algorithm  results  for  a up  to  2.5  where  stall  took  place.  The  pressure 
distributions  over  the  airfoils  and  wake  centerlines  in  the  trailing  edge 
regions  are  shown  in  Fig.  25.  As  the  trailing  edge  angle,  a increases,  the 
wake  effect  is  anticipated  further  forward  on  the  airfoil  surface  in  much  the 
same  manner  as  in  the  compression  corner  case.  A detailed  comparison  of 
Figs.  22  and  25  indicates  that  the  interaction  effect  is  slightly  less  severe 
for  the  trailing  edge  problem  then  it  is  for  the  compression  corner.  The 
airfoil  surface  shear  distributions  for  the  region  ahead  of  the  trailing  edge 
are  given  in  Fig.  26  for  the  full  range  of  a studied  here.  Comparison  with 
Daniels  (Ref.  11-14)  flat  plate  results  are  seen  to  be  very  good  even  in  the 
immediate  trailing  edge  region  where  the  rapid  acceleration  in  the  wake  region 
has  caused  a thinning  of  the  viscous  layer  and  subsequent  increased  shear.  As 
the  trailing  edge  angle  was  increased,  the  initial  acceleration  effect  of  the 
wake  region  (producing  increased  surface  shear)  changed  to  a deacceleration 
effect  (producing  decreased  surface  shear).  The  appearance  of  negative  shear 
for  a * 2.5  signals  the  occurrence  of  reverse  flow  (i.e.,  a stalled  trailing 
edge).  The  accompanying  wake  centerline  velocity  distributions  for  the  same 
range  of  a are  given  in  Fig.  27.  Here  comparison  with  Daniels  (Ref. 11-14) 
flat  plate  results  show  a slight  loss  in  accuracy  apparently  due  to  the 
extremely  rapid  flow  accelerations  experienced  at  the  trailing  edge  (Daniels 
shows  that  — * • as  X ♦ 0+).  Fortunately,  as  the  trailing  edge 
angle  was  ift£reased,  this  acceleration  decreased  until  it  finally  reversed, 

. indicating  that  a recirculating  stalled-flow  bubble  was  formed  for  a ■ 2.5. 

It  was  found  that  the  computing  time  required  for  the  second  order  algo- 
rithm to  achieve  a converged  solution  increased  as  the  trailing  edge  angle 
increased.  Whereas  the  solutions  for  a * 0 were  achieved  in  1.5  minutes  of 
CPU  time,  approximatel y 9 minutes  were  required  for  the  a ■ 2.5  case.  Since 
the  first  order  algorithm  was  observed  to  achieve  convergence  aignif icantly 
faster  than  its  second  order  counterpart,  solutions  were  also  obtained  here 
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over  a wide  range  of  trailing  edge  angle  using  this  algorithm.  The  resulting 
pressure  distributions  are  shown  in  Fig.  28  for  a up  to  3.  These  results  are 
essentially  the  same  up  to  a = 2 as  those  given  in  Fig.  25  for  the  second 
order  algorithm  but  were  obtained  here  in  approximately  one  third  of  the  com- 
puting time.  For  a = 2,  the  difference  in  the  results  was  found  to  be  small 
and  to  diminish  as  the  step  size  was  reduced  in  the  first  order  algorithm. 
Thus  this  algorithm  provides  a reliable  and  reasonably  efficient  method  for 
solving  stalled  trailing  edge  flows  and  should  be  useful  in  the  extension  of 
the  present  approach  to  other  flow  regimes  (i.e.,  subsonic  flows,  loaded 
airfoils,  turbulent  flows). 
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CONCLUDING  REMARKS 


Two  major  areas  have  been  investigated  under  the  present  effort.  The 
first  concerns  the  application  of  a numerical  procedure  for  solving  the  com- 
pressible Navier-Stokes  equations  to  cascade  blade  passage  flows,  and  the 
6econd  deals  with  the  construction  and  numerical  solution  of  an  interacting 
boundary  layer  model  for  airfoil  trailing  edge  flows. 

Under  the  first  part  of  the  present  effort  an  existing  computer  code, 
based  on  an  alternating  direction,  implicit  time  marching  solution  of  the 
general  coordinate  representation  of  the  Navier-Stokes  equations,  has  been 
modified  to  be  applicable  to  the  blade  passage  problem.  In  particular,  the 
general  coefficient,  matrix  inversion,  boundary  condition,  and  spatial  dif- 
ferentiation subroutines  of  the  existing  code  were  modified  to  incorporate 
cascade  blade-to-b lade  periodicity  conditions.  The  latter  consist  of  obtain- 
ing an  explicit  solution  of  the  governing  equations  at  one  periodic  boundary 
and  specifying  periodicity  of  density  and  velocity  at  the  other  periodic 
boundary . 

Implicit  time-marching  Navier-Stokes  solutions  have  been  computed  for 
flows  past  a circular  cylinder  at  Reynolds  numbers  of  forty  and  eighty.  The 
convergence  rate  for  this  code  was  found  to  be  slow.  In  particular,  160  time 
steps  for  the  Re  = 40  case  and  200  time  steps  for  Re  * 80,  with  At  varying 

_ *) 

between  0.15  and  0.5,  were  required  to  achieve  a convergence  level  of  2 x 10 
(c.f.  Eq.  (1-34).  In  addition,  for  the  Re  * 40  case  pressure  and  velocity 
distributions  after  160  time  steps  are  not  in  good  agreement  with  previous 
calculations  or  experimental  data.  Time  history  plots  of  the  current  calcu- 
lations indicate  that  continued  iteration  would  not  improve  the  result  compar- 
isons. At  present  the  reasons  for  the  lack  of  agreement  between  present  and 
previous  calculations  are  not  apparent,  and  further  work  (e.g.,  step  size 
studies  to  interrogate  truncation  errors,  investigation  of  different  coordinate 
distributions,  etc.)  seems  warranted  to  resolve  this  is6ue. 

Under  the  second  part  of  the  present  effort  an  interacting  boundary 
layer  model,  incorporating  the  asymptotic  triple  deck  concept,  has  been 
formulated  for  airfoil  trailing  edge  flows.  Attention  is  focused  on  the 
flow  past  symmetric,  wedge-shaped  trailing  edges,  with  trailing  edge  angles 
large  enough  to  cause  flow  separation  from  the  airfoil  surface  and  re-attach- 
ment on  the  wake  slip  line.  The  analytical  model  is  taken  directly  from 
that  already  proved  to  be  applicable  for  flow  into  a compression  corner. 

Numerical  solutions  to  the  triple  deck  version  of  this  model  have  been 
obtained  for  the  laminar  supersonic  case.  Results  determined  for  a 0*  trail- 
ing edge  angle  (flat  plate)  were  found  to  be  in  good  agreement  with  those 
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obtained  earlier  by  P.  G.  Daniels.  Present  calculations  have  been  extended 
over  a range  of  trailing  edge  angles  until  a separated  flow  solution 
("stalled"  trailing  edge)  was  achieved.  Solutions  have  been  obtained  with 
two  different  finite  difference  formulations.  One  is  essentially  second 
order  accurate  while  the  other  is  first  order  accurate  but  provides  faster 
solutions.  Accuracy  studies  indicate  that  the  two  formulations  are  producing 
consistent  results  for  the  stalled  trailing  edge  case.  Although  both  attached 
and  separated  supersonic  trailing  edge  flows  have  been  sucessfully  calculated, 
there  are  areas  for  possible  improvements  in  the  present  numerical  algorithms, 
in  particular,  in  treating  flow  discontinuities  at  the  trailing  edge,  if 
experimental  comparisons  indicate  that  such  improvements  are  warranted. 

Future  work  should  also  be  directed  at  the  extension  of  the  current  approach 
to  turbulent  flows,  lifting  configurations  and  full  cascade  flow  fields. 
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Figure  8-Wake  Centerline  Velocity  Distribution  Behind  a Circular  Cylinder 
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Figure  17  Separation  Shear  Levels  for  High  Reynolds  Number  Flows 


78-08-97-1 


Figure  18  Separation  Pressure  Levels  for  High  Reynolds  Number  Flows 
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Figure  24  Accurac,  Details  in  Trailing  Edge  Region 
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Figure  26  Supersonic  Trailing  Edge  Shear  Stress  Distributions 
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